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Niunerical  Gclutions  are  obtained  to  the  sixtb  order  systooi  of 
disturbance  equations  describing  the  hydrodynaraic  stability'-  of  a 
laminar,  heated,  flat  plate  water  boundary  layers  including  all  mean 
and  disturbance  property  variations.   The  results  are  compared  v?ith 
those  calculated  using  the  fourth  order  system  of  Wazzan,  Okamura, 
and  Smith  [4]  (which  assumes  only  rriean  viscosity  variation). 

Both  systems  predict  sip;niflcant  boundary  layer  stabili2;ation 
(increased  minimum  critical  Reynolds  Number,  decreased  disturbance 
amplification  rates,  etc.)  with  moderate  heating,  but  display  a 
maximun  and  subsequent  decrease  as  the  plate  surface- to-free  stream 
temperature  difference  is  further  increased.   Over  thn  normal  liquid 
range  of  water,  the  results  of  the  sixth  order  calculations  shov;  only 
a  slight  enhancement  of  stability  over  those  obtained  from  the  fourth 
order  systera.  Apparent  gains  in  stability  predicted  V7hen  only  dis- 
turbance viscosity  terms  are  included  in  the  calculations  are  offset 
by  the  further  inclusion  of  denr.ity  fluctuation  terms. 

The  insensi tivity  of  stability  characteristics  to  the  addi- 
tional consideration  of  thermal  disturbances  is  attributed  to  the 
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very  lind.ted   region  in  v;hich  such   disturbances   are   important   for 
fluids   of   large  Prandtl  nuiuber  such  as  water. 
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CHAPTER  I 
INTRODUCTION 

To  gain  some  Insight  into  the  mechanism  of  laininar-turbulent 
boundary  layer  transition j  considerable  experinental  and  theoretical 
effort  has  been  expended  in  studying  the  bydrodynamic  stability  of 
a  laminar  boundary  layer.   This  work  has  revealed  that  alterations 
to  the  iTiean  flow  by  pressure  gradients,  boundary  roughness  and 
compliancy,  suction  or  blowing,  heating  and  the  like  can  cause 
significant  changes  in  stability  characteristics  of  the  flow.   The 
present  study  focuses  on  the  Influence  of  heating  on  the  stability 
of  a  liquid  boundary  layer,  and  is  motivated  by  the  favorable  effect 
that  surface  heating  may  have  on  drag  reduction  of  fitibtrio.rged  under- 
water vehicles. 

An  indication  of  the  qualitative  effect  of  surface  heating 
can  be  obtained  by  utilizing  the  inviscid  point-of-inflection 
criterion  (for  constant  density  flows  velocity  profiles  having  a 
point  of  inflection  are  unstable).   At  the  surface  of  a  flat  plate, 
the  boundary  layer  equations  for  a  fluid  flow  with  a  nonconstant 
viscosity  reduce  to 


del  ]    _    _    1      ( djJ\     f  du\ 


If  the  plate  is  heated  (such  that  T  >T^) ,  then  (dT/dy)  ^<0.   Because 
the  viscosity  of  a  liquid  decreases  v;ith  increasing  temperature 

(dp/dT  <  0),  then  (.d\i/dy)      >   0.   Conversely,  for  a  cooled  plate 

w 

(T  <  T  ),  (dp/dy)   <  0.   As  both  y   and  (du/dv)   are  positive  for 

W     '^^  W  W  '   w 

>       2    2   < 
either  case,  it  follows  directly  that  for  T    T  ,  (d  u/dy'")    0. 

2    2 
Since  the  curvature,  (d  u/dy  )j  is  vanishingly  small  but  negative 

as  y  ■*•  ",  it  follows  that  if  its  value  at  the  v/all  is  positive,  there 

must  be  at  least  one  point  of  inflexion  in  the  boundary  layer  (i.e., 

some  point  at  which  the  curvature  is  zero) .   Thus  it  is  anticipated 

that  heating  should  stabilize  and  cooling  destabilize  the  liquid 

laminar  boundary  layer.   This  conjecture  has  been  verified 

theoretically  by  a  number  of  authors.   (Note  that  it  is  just  opposite 

to  the  conclusion  reached  for  gases  where  viscosity  increases  with 

temperature) . 

DiPrina  and  Dunn  [l]  cite  unpublished  numerical  results  of 

Mcintosh  which  indicate  that  for  the  t%70-dimensional  Tollmicn- 

Schlichting  type  of  instability  in  a  heated  water  boundary  layer, 

the  minir.iuni  critical  Reynolds  number.  Re  .     .   ,  is  increased 

■'  min.crit. 

tenfold  over  the  isothermal  case  for  a  50  C  temperature  difference 
(T  =  10°C). 

^   CO 

Hauptmann  [2]  utilixcs  intep.ral  methods  and  a  perturbation 
procedure  to  correlate  variation  in  a  velocity  profile  shape  factor, 
arising  from  small  changes  in  a  tcm.perature  dependent  viscosity, 
with  the  stability  characteristics  for  a  class  of  single  parameter 
velocity  profiles.   His  results  also  predict  very  strong  stabilization 


in  water  for  relatively  smsll  heating,  and  are  quite  compatible 
vith  those  of  Mclntosli  cited  above. 

The  most  extensive  study  of  the  stability  of  a  heated  and 
cooled  v^ater  boundary  layer  with  a  Falkner-Skan  flow  has  been  con- 
ducted by  Wazsan,  Okaiiiura,  and  Smith  [3,4,5]  and  Wazzan,  Keltner, 
Okamura,  and  Smith  [6],   For  th.e  disturbance  equation,  assuming 
explicit  fluid  property  variation  in  onl^^  the  mean  viscosity  and 
no  fluctuating  temperature  field,  they  obtain  for  all  their 
analyses  a  modified  Orr-Sommerfeld  equation 

+ya"(0"^c^^(/j)l 

where  '  here  denotes  differentiation  with  respect  to  y/o.   In  con- 
trast, the  variable  raean-flow  coefficients  are  obtained  by  solving 
the  coupled  mean  momentum  and  energy  boundary  layer  equations  in 
which  all  fluid  properties  are  assumed  to  v^ry  [7],   In  general, 
their  results  show  that  while  cooling  the  wall  does  destabilize 
tlie  flow,  moderate  heating  very  strongly  stabilizes  it.   Hov/ever, 

as  heating  is  increased.  Re  ,     ,^   reaches  a  maximum  witli  T  and 

min  ,  cr:i  t .  w 

tlien  decreases.   They  explain  this  behavior  as  folloxvs  [s]: 


VThile  the  results  of  reference  [.3]  are  qualitatively  correct, 
they  arc  quantitatively  inaccurate  due  to  an  error  in  the  numerical 
integration  of  equation  (1.2)  and  subsequently  have  been  corrected 
[4,5].   The  first  reference  still  provides  the  best  explanation  of 
the  numerical  procedtirc  used. 


"As  T  is  increased  above  T  ,  the  velocity  profile  becomes 
more  stable,  whereas  the  variable  viscosity  tern.,  p  (  r\) ,    in  the 
modified  0-S  equation  tends  to  destabilize  the  flow.   However,  at 
moderated  rates  of  heating,  the  effects  of  the  nean  velocity  pro- 
files u  (ri)  ard  u"  (n)  upon  Re  .     .    are  the  dominant  ones.   As 

min.crit. 

the  rate  of  heating  increases,  the  negative  effect  of  y  (n)  upon 

Re  ,     .    increases  whereas  the  effect  of  heatinj;  upon  u  (n)  and 
rain.crxt. 

u"  (ii)  and,  in  turn,  their  respective  stabilizing  effects  on 

Re  ,     ,^  beein  to  level  off.   Re  .     .^   then  reaches  a  maximum 
rain.crit.   °  mm.crit. 

at  (T  )   ,   where  the  effects  of  y  (ti)  and  u  (n)    and  u"  (n)  upon 

Re  ,     ,   balance  out.  As  T  is  increased  beyond  (T  )   .   ,  the 
mm.crit.  w  •      v;  crit. 

destabilizing  effect  of  y  (il)?  becomes  dominant  and  Re  ,    _, 

decreases  upon  increasing  T  beyond  (T  )     ". 

Similar  results  and  conclusions  are  obtained  by  Potter  and 

Graber  [s]  for  the  stability  of  Plane  Poiseuille  water  flow  with 

heat  transfer.   Using  for  their  analysis  the  same  equation  as 

Wazzan,  et.al.  and  neglecting  all  property  variation  otlier  than 

viscosity  in  the  solution  of  the  mean  equations,  they  found  that 

omission  of  the  mean  viscosity  gradient  terms  from  equation  (1.2) 

indicated  increased  stability  with  heating.   In  contrast,  inclusion 

of  these  terms  show  destabilization  with  increased  heating  (i.e., 

a  decrease  of  Re  .     ..    )• 
min.crit. 

Tlie  motivation  for  the  present  investigation  stems  from  an 
attempt  to  more  fully  explain  the  stability  anomalies  found  by 
Wazzan,  Okanura,  and  Smith  through  a  more  complete  reformulation 


of  the  problem.   Specifically,  it  is  unclear  X'/hy  those  authors 
solve  the  coupled  mean  flow  system  of  nioraentum  and  energy  equations 
to  provide  the  variable  coefficients  for  the  disturbance  equation 
(1.2),  and  yet  neglect  a  disturbance  energy  equation  coupled  with 
the  momentum  equation  through  fluctuating  fluid  property  variation. 
As  an  a  priori  estimate  of  the  relative  magnitude  of  all  disturbance 
quantities  and  their  derivatives  is  impossible  throughout  the 
entire  boundary  layer,  the  omission  of  temperature  fluctuations  in 
previously  cited  investigations  is  seemingly  arbitrary  and  without 
justification  in  a  situation  involving  heat  transfer.   Accordingly, 
the  complete  linearized,  parallel  flow  disturbance  equations  (ex- 
cluding only  negligible  expansion  work  and  viscous  dissipation) 
eventually  will  be  formulated  and  solved. 

A  simple  illustration  substantiates  the  preceeding  conten- 
tion.  Assu-ning  for  the  moment  that  the  only  property  variation  of 
importance  is  that  of  viscosity,  it  will  be  shown  that  the  stability 
equations  can  be  written  more  completely  as 


(1.3) 


coupled  v;ith  a  disturbance  energy  equation 


Lt(u-c)  +  hV  =     V  p 

-^      '^     o,'  Re  I  r 


(1.4) 


through  the  relation  betv/een  the.  viscosity  fluctuations,  m,  and 
the  temperature  fluctuations  t .   Assuming  further  that  disturbance 
quantities,  m  and  <f:,    as  well  as  their  derivatives,  m"  and  ({>",  are 

each  of  comparable  magnitude  at  some  point  in  the  boundary  layer 

t  

and  realizing  that  u  is  of  the  same  magnitude  as  y",  then  terms 

arising  due  to  the  presence  of  the  disturbance  viscosity  such  as 

in"u  and  'X'^mu  are  likewise  of  comparable  magnitude  to  terms 

retained  by  Wazzan,  ct.al.such  as  (J)"vi"  and  o<'"(*,-,j"  ,  and  so  should 

not  he  discounted  in  the  analysis.   Since  V.'azzan,  et.al.  attribute 

much  of  the  stabilisation  to  the  mean  viscosity  gradient  terms  [s], 

it  would  seem  as  if  omission  of  these  terms  is  "not  v/ell  justified" 

and  it  is  not  petTnlssible  categorically  "to  assume  that  T  =■•  0 

throughout  the  boundary  layer"  [3]. 

Thus  the  focus  of  the  investigation  is  upon  establishing  the 

difference  between  the  fourth  order  system  of  Uazzan,  et.al., 

equation  (1.2),  and  a  more  complete  version  of  the  sixth  order 

system  represented  by  equations  (1.3)  and  (1.4)  in  v;hich  all 

property  variation,  mean  and  disturbance  (acceptable  within  the 

locally  parallel  flow  and  boundary  layer  assumptions),  is 

considered. 


A  complete  discussiori  of  the  assumptions  made  in  formulating 
the  problem  considered  in  this  investigation  is  given  in  Chapter 
II.   The  assumed  fluid  property-teniperature  relationships,  niunerical 
techniques,  and  assessir.ent  of  the  numerical  errors  incurred  in  the 
process  are  discussed  in  Chapter  III.   The  results  of  this  analysis 
are  presented  in  Chapter  IV,  and  finally  the  conclusions  dra^cn  from 
these  results  are  discussed  in  Chapter  V.   The  Appendices  are  in- 
cluded to  present  relevant  but  nonessential  naterial  that  would 
otherwise  obscure  the  primary  thrust  of  this  investigation  in  the 
text. 


CHAPTER  II 
FORMULATION  OF  THE  STABILITY  PROBLEM 

This  .investigation  examines  the  spatial  stability  to  small- 
scale  disturbances  of  the  steady,  laininax^,  boundary  Jaycr  flow  of 
a  viscous,  heat-conducting  liquid.   Since  formulation  of  the 
problem  can  easily  be  nadc  with  far  greatei"  generality  than  is 
eventually  required,  such  a  procedure  is  follot.'ed  to  demonstrate 
that  v/ithiu  the  assuraptions  to  be  made,  the  additional  coraplexities 
fail  to  significantly  alter  the  resulting  stability  equations. 

Thus  while  the  problem  is  posed  for  general  Falkner-Skan  flows 

M 
(u  =  Cx  )  V7ith  constant  cross  flows  (v?  =  constant)  and  ax"bitrary 
e  e 

disturbance  orientation,  solutions  discussed  in  Chapters  VI  and  V 
are  specifically  for  tuo-dimensional,  flat  plate  boundary  layers 
with  tV70-dimsnsional  disturbances.   Tiie  Falkner-Skan  flovjs  are 
often  referred  to  as  "v,'edge  flovvs"  since  for  M  >  0  they  correspond 
to  a  uniform  How  over  a  yavjed,  rigid,  infinite  v.'edge  v.'ith  included 
angle  "rrrj —  [9,  p  lAl]  at  zero  angle  of  attack  in  the  neighbor- 
hood of  the  stagnation  point. 

2.1  BasjLc  Assuirptions  and  Goyerning  Equations 

To  describe  this  flow  field  mathematically,  certain  basic 
assumptions  are  made  concei:ning  the  nature  of  the  f].uid  involved. 

(1)  'Iho  characteristic  scales  for  mean  motion  and 


disturbances  to  it  must  be  signif icaatly  greater  than  those 

associateid  vrith  molecular  motions  in  the  fluid  (continuum 

hypothesis).   This  pcnrdts  ascription  of  meaningful  values  to 

fluid  and  flow  parameters  "at  a  point".   Further,  for  disturbance 

periods  (2i;/w)  much  greater  than  the  molecular  tine  scale,  say 

the  fluid's  structural  relaxation  time  (T  )  (i.e.  wT  <<  1),  the 

s  s      ' 

flow  can  be  considered  to  be  in  statistical  equilibrium  so  that 
molecular  transport  effects  can  be  represented  by  transport  co- 
efficients (v!.  k,  etc.). 

(2)  Once  the  continuum  hypothesis  is  established,  to  define 
the  stress  relations  bct^-7cen  fluid  elements,  it  is  assumed  that 

the  fluid  is  Newtonian  (linear  stress  variation  with  rate  of  strain) 
and  isotropic  (no  intrinsically  preferred  elemental  spatial 
orientation) . 

These  two  groups  of  assumptions  form  the  foundation  for  the 
Navier-Stokes  equations  (2.1-2.3).  Since  there  is  no  experimental 
evidence  to  suggest  their  inadequacy  under  ordinary  conditions  , 

they  will  be  used  without  reservation  in  this  analysis. 

(3)  Ass'oming  the  fluid  to  be  a  liquid,  the  form  for  the 
state  equations  is  defined.   The  relative  pressure  independence  of 
all  fluid  properties  permits  their  expression  explicitly  as 
functions  of  ter.peraturc  alone. 
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Thus,    (dinc-asionalljO 


Continuity: 


D 


Dt 


^  ^  ^^K.      ■=     O 


(2.1) 


Moment  urn: 


Energy : 


State: 


V  Dt 


C 


6.  4-  M, 


c)\ 


TJ 


D" 


D 


P.  C   , 


/v 


-Funct 


ion 


(2.2) 


(2.3) 


(2.  A) 


where      :e-  -  material  derivative  =   —  +  L^,-  .— 


at    J  ^^ 


Cr  =  stress  tensor  =  -  p|.  S^c,'  4-  2/x€,-;  4-  ?\^^^($ 


^J 


£,•,■  =  rate-of -strain  tensor  =  -r  I    ^=^'  +  ^^^ 


[3  =  body  force 


^.  4-  All  ^ 


S  ~  volumetric  expansion  coefficient  =  -  —  (  -J- 

P~Pt      -   "mechanical"  ~  thermodynamic  pressure 
difference  =  -  (>^|/v>)£^^ 


f 


This  last  equation  relating  the  mechnnical  pressure  (negative 

of  the  mean  normal  stress,  -  ~  cr,-,-)  arid  tlicrmodynamic  pressure  is 

valid  because  only  the  near  past  is  considered  (or  wT  <<  1  as 

s 

previously   discussed    [lO]    ).      l\ote   that   equations    (2.2)    reduce   to 


the   Navier-Stokcs    equations    if  — 


(Stolceij'    liypothesis, 
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valid  only  for  an  ideal  monatomic  gas)  or  t,,,  =  0. 

For  liquid  in  a  heated  boundary  layer,  there  is  no  justifi- 
cation for  assuming  either  statement  is  valid  in  general. 

Liebermann  [11 J  '^^s  shown  that  for  a  number  of  liquids ^  —  is 

y 

positive,  considerably  greater  than  vmity  in  some  cases,  and  widely 
different  for  different  liquids.   Furthermore,  for  low  frequencies 
(wT  <<  1),  NarasirrJiara  [l2]  notes  that  the  viscosity  ratio  is 
independent  of  disturbance  frequency,  and  almost  independent  of 
temperature  [l3].   (For  water,  Pinkerton  [1A]  has  shovm  a  — 
variation  of  about  16%  over  the  range  0-60  C  as  compared  with  a 
y  variation  of  118%  over  the  same  range  -  see  Figures  B.2  and  B.6). 
Alternatively,  for  flows  v^ith  heat  transfer  (particularly  cases 
with  large  thermal  gradients),  when  the  density  is  identified  as  a 
function  of  temperature,  it  will  not  be  constant  in  that  flow 
field.   In  such  a  situation,  then  e   will  vanish  only  at  a  rigid 
surface  for  steady  flows.   Thus  if  relative] y  large  amplitude 

density  fluctuations  were  expected  in  a  liquid  for  which  was 

large  (e.g.,  for  benzene,  this  ratio  is  0(10  )  [ll],  then  there 
could  be  a  large  disparity  between  mechanical  and  thermodynamic 
pressures.   It  will  be  shovm  in  Appendix  A  and  Section  2.3.1  that 

for  values  of  of  order  unity,  the  stability  problem  can  be 

solved  independently  of  this  quantity. 

l«rhen  the  flov;  is  characterized  by  a  single  time  scale,  namely 
that  associated  with  the  frequency  of  the  small,  rapid  perturbations 
on  the  basic  flow,  cacli  of  tlic  flow  and  transport  parameters  can 
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be  constructed  as  the  sum  of  a  time-averaged  mean  and  a  fluctuating 
coir.poneri.t. 


Q.(VO=  Q.(!,)^Q'^,,i 


)  (2.6) 


Equations  describing  the  behavior  of  these  components  are  developed 
in  sections  2.2  and  2.3.  respectively. 

2 . 2  The  Mean  Flow 

The  ir.ean  equations  for  boundary  layer  flow  over  a  v/edge  can 
be  deduced  from  equations  (2.1-2.3)  (see  Appendix  A)  after  choosing 
the  orthogonal  coordinate  system  shown  iii  figure  2.1  and  making  the 
following  simplifying  assumptions: 

(1)  The  basic  f3ow  is  steady  but  not  necessarily  spatially 
homogeneous  (cqn.  2.6),  so  the  mean  equations  evolve  from  a  time- 
average  of  the  instantaneous  equations. 

(2)  In  the  momentum  equations  for  motion  parallel  to  the 
surface,  viscous  and  inertial  forces  are  of  the  same  order  of 
magnitude.   Using  the  scaling  resulting  from  this  requirement  and 
eliminating  terms  of  0(Re-   )  or  smaller,  boundary  layer  equations 
are  obtained. 

(3)  The  wedge  is  "infinite"  in  extent  in  the  z   direction,  so 
all  jnean  flovj  variation  is  independent  of  the  z  coordinate. 

(A)  Gravity  is  the  only  body  force  present,  acting  in  a 
plane  normal  to  the  z  axis  at  an  angle  x  to  the  y  axis. 
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Additional  siTr.plifications  leading  to  exclusion  of  specific 
effects  on  the  flov?  (e.g.,  dissipation,  compression  energy,  buoyancy) 
will  be  made  after  examining  relevant  nondimensional  coefficients. 

2.2.1  Mean  Flov.''  Equations 

Using  the  assuir.ptlons  discussed  in  sectioxi  2.2,  equations 
(2.1-2.3)  reduce  to  the  (time-averaged)  mean  boundary  layer 
equations  (din;»usionally) 


Continuity:  ^4^  4-  ^_eJ:t  =  Q  (2.7) 

by  dy 

Momentum:     p  f  ij   ^^  +   V   ^^  ]  ^    -  P  °\   ^"-^^  -   '^f*         (2.8) 


+    ^     l/A 


ay  - 


V  ^   ^x      by  /        c>yy     ^)y/ 


Energy:    7.  C  f  ^  --  4-  v  ii  )  t:  (b  1  ^  ^Lft  (2.11) 


y/   V  oy/  J  by  ^      dy 

Before  solving,  several  modifications  are  made  to  the  x 
momentum  equation.   Evaluation  of  this  equation  at  the  boundary 
layer  edge  reveals  the  nature  of  the  mean  pressure  variation  im- 
pressed on  the  boundary  layer  (dimensionally) 
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f  =  ^»  --=  f  .  e. 


Sin 


-y 


(2.12) 


Further,  due  to  the  state  equation  (2.4),  the  density  can  be  ex- 
panded in  a  Taylor  series  in  temperature,  so  dimensionally 


'    1  et 


+ 


or  nondimensionally,  using  the  volumetric  expansion  coefficient, 
6,  defined  in  equations  (2.5), 


l-p( 


^^i)=  (^^X.  (T-  Oj  1 


+ 


Ct-1) 


z 


+-  ftccTo 


For  small  values  of  B  T  and   «>  -  l   (such  as  is  the  case  for 

T 

CO 

V7ater  with  noderate  temperature  differences)  ,  this  can  x«>ll  be 
approximated  by 


1  -   P  (t) 


((t)-  ^^  U(T-1) 


(2.13) 


Thus,  with  the  modifications  indicated  by  equations  (2.12) 
and  (2.13),  the  x  moraentum  equation  finally  becomes  (dimensionally) 


^ 


00    ^e 


o/x 


+  Poo  ^ooq  sinX  (1  -"u) 


+-     ^-     ^u    ^^^ 


dy  y     dy 


(2.8a) 
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Since  all  fluid  properties  are  asEumed  to  be.  variable, 
particularly  density,  a  modified  Uowarth-Dorotnitsyn  transformation 
is  suggested  [7],  reducing  the  boundary  layer  equations  to  a 
nearly  Incompressible  form.   The  required  independent  variables 
are 


X-  X 


> 


rfi~r]  i'y 


y  (2.14) 


Satisfaction  of  the  continuity  equation  is  accomplished  by  intro- 
ducing  a  dirr^ensional  straam  function  ^     such  that 

p  u     -     i- 

e^v'  =  -  ^'^^ 

and  a  nondiDiensional  stream  function  F  so  that 


where  the  two  are  related  by 

Defining  the  displacement  thickness 


^/T'  (2.15) 


^'-~      r  ()-   ^u)cJy^ 


(2.16) 
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or  in  transformed  variables 


-  U 


J 


7 


(2.16a) 


then 


where 


^ 
h* 


(2.17) 


y 


I 


uU 


1 


(2.18) 


a  constant  for  specified  v/all  and  free  stream  temperature  and  wedge 
angles. 

Additionally,  defining 

w*  =  Wg  G  (n) 
then  equations  (2. 7. 2. 8a, 2. 10,  and  2.11)  transform  to 

D    -      /.  .      .  r 

-    PrE  j   P/l      F 


Re  ^ 


/  n' 


\  V       2 


^  /!iM  G'"U  1/-^-" 


t/e/ 


Gr 


'w-'cv      I?; 


•Sin  \ 


(2.19) 


Vc 

r 

i/e 

^rJ 

CO 

c^. 

/t/» 

fie  9,  C"^v.-T^)^ 

'    V      ,,2^3 
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whera         Re  =  Reynolds  Number  = 

Gr  =  Grashof  Number  = 

E  =  Eckcrt  Number  = 

Pr  =  Prandtl  Number  =   

M  =  — -  ^      =  constant  since  u  is  of  tne 

Ue   dx  e 

M 
foni  u  =  X  for  the  wedee  flov;. 
e 

At  this  point  an  assessment  is  made  of  the  various  force  and 
energy  terms  of  the  above  equations  to  determine  their  influence 
in  establishing  the  mean  velocity  and  temperature  boundary  laj'er 

profiles,  vrith  the  following  conclusions. 

Gr 

(1)  Buoyancy  effects  can  be  neglected  if  -- —  sin  x  "^'^  !• 

Ke 

This  is  easily  accomplished  analytically  by  requiring 

(a)  the  surface  normal  be  parallel  with  the  gravity 
vector  so  that  x  =  0,  or 

(b)  for  large  x  angles  (i.e.  sin  x  ==  0(1))  ,  u  be 

Gr 
sufficiently  large  that  r— -  <<  1,  since  this 

Re 

2 
parametric  group  varies  inversely  vrith  u  .   Thus, 

except  for  surface  orientations  nearly  parallel  to 

the  gravity  vector  and  low  speeds,  free  convection 

can  be  neglected. 

(2)  Disf;ipation  and  expansion  (or  pressure)  vjork  are 
negligible,  requiring  that  Pr  E  <<  1.   (An  estirate  of  these  terms 
for  a  water  flow  indicate  that  they  can  never  be  significant  while 
the  boundary  layer  approximation  remains  valid) , 
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As  a  consequence  of  these  assumptions,  the  boundary  layer 
equations  (2.19)  are  explicitly  independent  of  x.   The  original 
set  of  partial  differential  equations  has  been  transformed  to  one 
of  ordinary  differential  equations,  so  the  velocity  and  temperature 
profiles  are  given  by  similarity  solutions. 

The  laean  flow  equations  can  be  ^^7ricte^  finally  as 


, F-y  .  M  (4  -  f'').(ma_')ff"=o 


B^G')'W^')FG'=0 


SRW')'  +  Prc(^)FI-{=0  (2.20) 


Note  that  the  z   moinentuin  equation  is  not  coupled  VTith  the 
other  two;  that  is,  F  and  H  are  determined  independently  of  G 
provided  p,  p,  and  k  arc  also  independent  of  G. 

2.2.2  Mean  Flow  Boundary  Conditions 

The  boundary  layer  equations  in  section  2.2.1  are 

supplemented  with  the  following  boundary  conditions: 

^       y 

aty  =   0  u=     w=0        (no   slip   at    the  v/all) 

V  =   0  (impermeable  v;all   surface) 


T*=  T 
w 


(constant  vzall   temperature)    .„    ^-^ 
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as  y  ->■  "    u  -HI  (x) 
^  e 


— )* 

e 


T  ->T 


Note  that  although  u  may  vary  with  x  due  to  a  nonzero 
v;edge  anple,  the  free  stream  temperature  (identical  with  that  at 
the  boundary  layer  edge)  is  assuvaed  constant. 

In  terms  of  the  transformed  variables,  these  conditions 
become 

atn=0      F'=G=0 
F  =  0 

"  "  ■*■  (2.21a) 

as  n  ■*  "•     F'  ->■  1 
G  -^  1 

H  -+•  0 

2 • 3  The  Disturbance  Flow 

The  disturbance  equations  are  obtained  as  the  difference  be- 
tween the  instantaneous  and  mean  equations  (see  Appendix  A) . 
Simplification  is  accomplished  by  the  following  considerations. 

(1)  Consistent  with  the  assumptions  of  section  2.2, 
"infinitesimal"  disturbance  amplitudes  permit  linearization  of  the 
disturbance  equations.   (Mote  that  all  term.s  of  this  order  vanished 
in  taking  the  mean).   This  of  course  eliminates  terms  with 
correlated  disturbance  quantities  (as  for  the  mean  equations)  and 
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reduces  the  disturbance  equations  to  a  set  of  linear,  partial 
differential  equations  v;ith  nonconstant  mean  flow  coefficients. 

(2)  The  mean  boundary  layer  flow  is  assumed  to  be  locally 
parallel,  so  that  for  a  specified  distance  from  the  leading  edge, 
X,  these  coefficients  of  the  disturbance  quantities  are  functions 
only  of  the  distance  normal  to  the  solid  surface,  y.   That  is, 
the  stability  of  a  velocity  profile  at  a  given  x  is  to  be  con- 
sidered independently  of  the  rest  of  the  boundary'  layer.   The 
requirements  necessary  to  make  such  an  assumption  are  elucidated 
in  Appendix  A  and  briefly  summarized  below.   ^-Jhen  the  amplitude  of 
all  disturbances  and  derivatives  of  disturbances  are  equal  (though 
not  necessarily  equal  to  each  other) ,  comparison  of  the  mean  flow 
coefficients  in  a  boundary  layer  sense  indicates  that  v  <<  u  , 

w  and  ■■■°-  -  <<  — ;-i— •  .   Furthermore,  if  disturbance  wave  lenc;ths 
3x       dy 

are  small  compared  v/ith  lengths  characterizing  changes  with  x  of 
mean  quantities,  then  these  quantities  can  be  considered  to  be 
slowly  varying  functions  of  x  compared  with  the  fluctuations. 

Again,  the  disturbance  equations  are  formulated  as  generally 
as  possible  within  the  fram.cvjork  of  assumptions  specified  above. 
Elimination  of  other  terms  by  examining  the  magnitudes  of  their 
dimensionless  coefficients  vjill  lead  to  additional  simplifications. 

2.3.1  Disturbance  Equations 


With  the  assumptions  discussed  in  section  2.3  the  disturbance 
equations  for  this  problem  become  (dimensionally) 
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Energy; 
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Notice  that  by  virtue  of  the  parallel  flow  assumptions,  the 
mean  flow  coefficients  of  all  disturbance  terms  are  functions  of 
y  only,  so  these  equations  are  noiJ  separable  in  terms  of  solutions 
of  the  form 

Noting  that  the  assuraed  disturbance  form  is  complex  (x-?ith  complex 
amplitude  function  q)  while  the  disturbance  equations  (2.22-2.26) 
are  not,  physical  quantities  can  later  be  recovered  by  taking  the 
real  part  of  the  conplex  solution.  Thus,  the  disturbances  are 
assumed  to  be  sinusoidal  waves  of  frequency  to  =  o<  c  propagating 
at  an  angle  0  to  the  x  axis.   For  this  spatial  analysis,  w  is 

•A 

assumed  real  and  d<  complex  to  admit  solutions  which  are  temporally 

A  A 

only  oscillatory  but  spatial  ly  amplified  (c<   <  0)  ,  neutral  (g(.  =  0)  , 


2A 


or  damped  (c<^  >  0)  , 

Such  a  formulation  vjxll  reduce  the  set  of  disturbance 
equations  froin  partial  differential  equations  to  ordinary  ones  in 


y.   Thus  introducing 
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v;here 
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and  defining  for  the  mean  and  disturbance  velocities  parallel  to 
the  surface,  respectively, 
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where  tan  Y  = 


and  6  is  the  angle  between  the  x  axis  and  the 


direction  of  v.'ave  propagation,  then  the  disturbance  equations 
become: 
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v/here  differentiation  is  again  v:lth  respect  to  n  to  permit  in- 
formation obtained  from  the  mean  equations  to  be  used  directly  in 
the  numerical  solution  of  the  disturbance  equations,  and  where 
all  fluid  property  variation  is  of  the  form 
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The  dimensionless  wave  number  and  Reynolds  number  are  defined  by 


o(     = 
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and 
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Note  that  the  composite  x-z  momentum  equation  is  obtained  by 
multiplyin;;  the  transformed  x  and  z  momcntura  equations  by  cosO  and 
sine,  and  then  adding.   Furthermore,  diJatory  expansion  energy 
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drops  out  of  thfi  energy  equation  with  li.nearization,  and  so  must 
be  considered  a  second  order  effect  for  this  problem. 

To  simplify  these  equations  further,  the  assumptions  made  in 
section  2.2.1  regarding  buoyancy,  dissipation  and  expansion  energy 
are  used  again,  with  some  interesting  iinpj.ications. 

(1)  Neglect  of  the  dissipation  energy  terms  (for  '^r^    -^"^    \      ) 
permits  the  energy  equation  to  be  written  as  a  function  of  J     and 
W  rather  than  f  j  h,  u,  and  u,  thereby  rcdiicing  the  order  of  the 
system  of  equations  by  tv.'o  from  eight  to  tix.   In  other  words,  v;hen 
viscous  dissipation  is  negligible  in  the  disturbance  energy  equation 
then  "the  stability  of  a  three-dimensional  boundary  layer  to  a 
plane-wave  disturbance  of  arbitrary  orientation  reduces  to  a 
two-dimensional  stability  problem  governed  by  the  boundary-layer 
profile  in  the  direction  of  wave  propagation  and  by  the  mean  tem- 
perature profile"[15]  .   Thus,  alterations  of  the  external  flo\<j   due 
to  changes  in  pressure  gradient  or  cross  flow  v/ould  be  manifest 

as  changes  to  the  variable  coefficients  in  the  disturbance  equations. 

(2)  Iftien  the  expansion  energy  is  negligible  (B^^T^E  <<1)  , 
variation  of  the  tliermodynnmic  pressure  appears  only  in  the  momentum 
equations.   Under  these  circumstances,  explicit  inclusion  of  the 
second  viscosity  coefficient,  X  ,  in  the  momentum  equations  can  be 
avoided  by  using  the  "mechanical"  rather  than  thermodynamic  pressure. 
Such  a  change  would  not  affect  the  mean  pressure,  for  using 
equation  (2.j)    and  the  "parallcl-f lov?"  and  linearity  assumptions, 

p  =  p  ;  however,  it  v.'ould  result  in  a  difference  between  the  two 
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fluctuation  pressures    P  ~  Pt  ~  I  ■^   3/  /^  by.       hy        dzj 
or,  in  terms  of  equations  (2.27) 

(w'-  c) 


O 


^-^."-^ 


f    P0  +    tr   ^-— ^;  (2.36) 


As  noted  before,  for  the  large  Reynolds  and  small  V7ave  numbers 
anticipated,  this  pressure  difference  could  only  be  significant 
for  large  values  of  X/y**. 

On  the  other  hand,  if  the  expansion  energy  is  not  negligible 
(say,  (B  T  )b~0(l)),  the  appearance  of  X  can  still  be  avoided  in 
the  momentum  equations.   ^Elimination  of  the  thermodynamic  pressure 
fluctuations  between  the  tvro  momentum  equations  (  2.33) 
yields  a  vorticity  equation  -  an  Orr-Sommerfeld  type  equation 
modified  by  the  inclusion  of  fluid  property  terms  -  in  which  terms 
containing  X  identically  cancel  [16].   However,  this  procedure 
would  require  evaluation  of  higher  order  temperature  derivatives 
of  the  fluid  properties.   Since  the  equations  will  ultimately  be 
reduced  to  a  system  of  first  order  equations  for  solutioi^,  such  a 
step  should  be  avoided  if  possible. 

In  short,  the  second  viscosity  coefficient  should  influence 
only  the  disturbance  thermodynamic  pressure  distribution.   For 
water  flow,  it  has  already  been  noted  that  E  must  be  very  much  less 
than  unity,  so  the  first  procedure  for  eliminating  X  is  utilized 
for  this  analysis. 

(3)  Buoyancy  effects  are  negligible  ( ^r^-  <<  1). 

Si'  Ke 
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Since  the  assumption  vas  made  in  solving  for  the  mean  equations 
that  buoj^ancv  terms  of  order  Gr/Re  <<  1,  it  is  certainly  reasonable 

to  assure  that  the  body  force  influence  v;ith  terms  nov.'  of  order 

_  2 
(Gr/«Re  )  is  even  less  significant.   If  buoyant  effects  were  to  be 

included,  their  contribution  would  appear  indirectly  through  modifi- 
cations to  the  mean  velocity  and  temperature  profiles.   Although 
it  is  knovm  that  for  flows  driven  primarily  by  body  forces,  heating 
from  above  stabilizes  v.'hile  heating  from  below  destabilizes  the 
flow  through  the  influence  of  thermal  instabilities,  this  effect 
will  not  be  considered  herein. 

Thus  with  the  assumptions  discussed  above,  the  disturbance 
equations  finally  become: 

Continuity;  cr  (W-c)  +    ^  "J^   +-  (z0)  -  Q  (2.37) 
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Energy: 
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Since  these  equations  are  not  in  a  form  suitable,  for 
numerical  Integration,  they  are  rearranged  and  recast  as  a  system 
of  six,  first-order  equations.   Denoting  the  dependent  variables 


as 
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the  six  equations  can  be  written  as 
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where  C.,  is  a  Gx6  coefficient  matrix  with  the  following  nonzero 
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Z'   is  obtained  fron  the  composite  x-z  raomcntum  ecuatlcn; 
Z'  ,  from  the  continuity  equation;  Z  '  ,  from  the  energy  equation; 
and  Z]    ,    from  the  y  raorientum  equation. 


2.3.2  Disturbance  Boundarv  Conditions 


The  disturbance  velocity  boundary  conditions  used  in  this 
analysis  are  easily  obtained  by  requiring  that  tliere  be  no  relative 
motion  between  the  fluid  and  the  solid  surface  at  this  interface 
(no  slip),  and  that  the  velocity'  normal  to  the  surface  must  also 
vanish  there  (impermeability).   In  notational  forai,  at  y  =  0, 

(J  =   w*-  0  (no  slip) 

\/   -    U  (impermeability) 

or  at  n  =  0,         ^^   0  -  (9  (2.43) 

The  thermal  boundary  condition  at  the  wall  requires  more 
careful  consideration  (see  also  discussion  by  Dunn  and  Lin  [17]). 
The  scaling  procedure  used  in  Appendix  A  indicates  that  for  a 
region  very  close  to  the  surface  (y/6  <<  1) ,  the  energy  equation 
reduces  to  the  one-dimensional,  unsteady  conduction  equation 
(dimensionally) 


The  inportant  point  to  note  is  chat  longitudinal  temperature 
fluctuations  are  negligible  in  this  region.   Thus,  solving  this 
equation  in  the  solid  and  requiring  that  disturbance  temperatures 
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and  heat  fluxes  be  continuous  at  the  liquid  -  solid  interface 
(y  =  0),  the  thermal  boundary  condition  becomes 
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and  the  subscript  "s" 
denotes  properties  of  the 
solid. 


It  is  simple  enough  analytically  to  require  that  the  thermal 
inertial  of  the  solid  be  sufficiently  large  that  thermal  fluctua- 
tions die  out  in  the  solid  very  close  to  the  surface.   Under  these 


conditions,  the  appropriate  boundary  condition  v;ould  be  t(o)  =  0 
boundary  condition  should  be  closely  scrutinized  when  applying 


(i.e. 


^.- 


^1  )•   Realistically,  however,  the  full 


this  analysis  to  a  particular  solid,  at  a  specified  frequency,  etc. 
For  this  analysis,  the  condition 


T(o^  =  0 


(2.A4) 


Is  chosen. 

Further  it  is  assumed  that  all  disturbances  are  bounded  at 


large  distances  from  the  surface 
y »-cD 
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2 . 4  Disturbance  Field  Energy  Balance 

Although  mathematically  the  stability  characteristics  of 
the  disturbance  flow  field  can  be  coinpletely  described  (xjithin 
the  fraFiework  of  specified  assumptions)  by  the  equations  formulated 
in  the  three  preceding  sections,  it  is  instriictive  also  to  develop 
the  theory  on  a  more  physical  foundation  as  well.   Insight  into  the 
physical  rr.echanisnis  at  v.'ork  in  the  boundary  layer  is  thus  sought  by 
examining  the  production,  dissipation  and  transfer  of  the  dis- 
turbance energies . 

The  appropriate  equations  are  obtained  by  first  multiplying 
each  of  the  x,  y,  and  z  disturbance  momentura  equations  (2.23-2.25) 
by  its  corresponding  disturbance  velocity,  u,  v,  and  w,  and  then 
suraming  the  three.   After  some  algebraic  manipulation,  a  meaningful 


equation  is  obtained  and  e>q)ressed  (dimensionally)  as: 
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As.'sutiiing  only   t\70-dimensional  mean  and   disturbance   flovjs, 
equation    (2.A6)    reduces   to: 
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The  time-averaged  version  of  equation  (2.47)  written  nondimens ion- 
ally  in  terms  of  disturbance  amplitude  functions  (see  equations 
(2,27)  is: 
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and  *  denotes  the  complex  conjugate. 


It  is  instructive  to  examine  each  term  at  this  point  to 
identify  the  physical  process  it  represents  and  assess  its  relative 
importance.   Thus,  the  indicated  terms  represent: 

(1)  the  rate  of  change  in  kinetic  energy  of  the  disturbance 
field.   Note  that  for  a  neutral  disturbance,  the  averaged  value 
is  zero  so  there  is  no  net  energy  change  in  the  boundary  layer. 
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(2)  the  energy  dravm  fror.i  the  mean  riiotion  by  the  working  of 
Reynolds  stresses  against  the  gradients  of  the  mean  velocity.   As 
this  is  the  only  source  of  energy  for  the  disti-irbance  field  of 
tVie  unheated  flow,  it  can  also  be  expected  to  be  significant  when 
there  is  heat  transfer.   Furthermore,  it  is  expected  to  be  greatest 
and  positive  in  the  region  roughly  between  the  wall  and  the 
critical  point,  signifying  energy  transfer  froin  the  nieaii  ^  the 
disturbance  flow,  but  small  and  quite  possibly  negative  from  there 
to  the  boundary  layer  edge. 

(3)  the  energy  extracted  from  the  mean  flow  due  to  the 
working  of  the  disturbance  viscous  stresses  against  the  nean 
velocity  gradients.   By  and  largej,  it  is  expected  to  have  the  same 
sign  as  term  (2)  and  so  should  represent  pri.mar;I  ly  a  source  of 
energy  for  the  disturbance  field.   It  should  be  raost  significant 
in  regions  where  temperature  fluctuations  are  largest  since  the 
viscosity  and  temperature  fluctuations  are  proportional. 

(4)  the  rate  at  which  the  pressure  fluctuations  do  work. 
Physically  this  energy  is  produced  by  two  different  mechanisms: 
(a)  the  periodic  expansion  and  dilation  of  each  fluid  element 
(caused  by  density  fluctuations)  in  the  presence  of  a  dicturbance 
pressure  field,  and  (b)  the  movcrent  of  this  same  fluid  element 
along  a  mean  density  gradient  also  in  the  presence  of  the  dis- 
turbance pressure.   Note  from  equatioTi  (2.48)  th.Tt  for  a  neutral 
disturbance,  the  contribution  from  the  first  mechanism  will  vanish 
at  both  the  v;all  and  the  critlcaJ  layer. 
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(5)  change  in  energy  resulting  from  a  nonzero  mean  viscosity 

gradient.   As  may  be  seen  from  equation  (2.48),  thin  tern  is  of 

order  §-     times  that  of  the  major  production  terms  (2). 

Re 

Since  the  general  o7rder  of  o(  decreases  and  Re  increases  with 

K 

moderate  heating,  it  is  expected  that  this  term  should  be  of 
negligible  significance. 

(6)  ire.chanical  energy  lost  from  the  disturbance  field  by 
viscous  friction.   This  loss  is  rcanifest  physically  as  a  heat 
source  which  in  turn  produces  minute  changes  in  the  disturbance 
temperature  (and  so  density)  level.   Since  the  term  is  negative 
definite,  it  will  alv7ays  represent  a  dissipation  and  t;hould  be  most 
important  near  the  vjall- 

(7)  energy  dissipation  from  thermal  fluctuations.  Specifi- 
cally, this  is  the  energy  lost  associated  v.'ith  the  dilation  of  the 
fluid  elements  due  to  their  expansion  and  contraction  V7ith  changes 
in  density  (caused  by  changes  in  temperature) .  This  effect  should 
be  small  in  liquid  boundary  layers  where  the  amplitude  of  dis- 
turbance density  fluctuations  is  relativejy  small,  even  with 
significant  mean  thermal  gradients. 

(8)  rate  of  transfer  of  pressure  and  dilatory  energy, 
vorticity,  and  disturbance  viscous  shear  stress.   If  the  disturbance 
field  is  neutral  (o<  =  0) ,  no  energy  is  transferred  parallel  to  the 
plate.   Further,  if  this  term  is  integrated  from  the  v;a]l  to  seme 
large  distance  (n  ->•  »)  ,  the  remaining  energy  contributions  will 
vanish  identically  due  to  the  homogeneous  boundary  conditions.   Thus, 
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the  term  acts  only  to  transfer  energy  froro  one  position  in  the 
boundary  layer  to  another,  and  for  a  neutral  disturbance  makes  no 
contribution  to  the  net  energy. 

Tv.'o  examples  are  nov;  examined  to  more  clearly  identify  the 
teriris  of  iiiterest  in  this  investigation.   First,  equation  (2.48) 
is  reduced  to  a  form  conpatible  v.'ith  the  assumptions  made  by 
Wazzans  et.al;  namely  in  the  disturbance  equations,  "p  =  constant 
(but  variable  in  the  mean  flov;)  ,  y  =  p  (y)  ,  and  t  =  0  (so  p  =p  =  0) 
Thus,  for  a  neutral  disturbance,  the  y-integrated,  nondimensional, 
energy  equation  is  composed  of  only  the  production  and  dissipation 
contributions: 


r n^ 
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=  - -'rR- ^'(^^^^4   -  ^-K/'(^^-fJ, 


i^ 
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(2.49) 


(Note  that  the  variable  p  still  appears  in  this  equation  since 
differentiation  is  with  respect  to  n ;  i.e.  A(\i''j^'-\        (V^   ^  X 

Secondly,  assuming  for  this  analysis  all  mean  and  disturbance 
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quantities  may  vary,  the  equivalent  equation  becomes 


(^) 


+   ^R 


(4) 


(5) 


^' 


-  Re  'Z('Tr)  -    45"  r/fpc^>l 


a) 


'A^  '^Vl 


"ZT" 


y^ 


'^-     ^      .(7)  ^ 


I-^ 
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t^^^(u~c) 


4-  4  ^^(p0) 


m  o(y 


IJ 


where   ^  (q)  is  here  slightly  simplified  to 


\^)-  -[g'°^':^'"''4  *  I'b 


',.  "■.  e'^ 


CHAPTER  III 
NUMERICAL  SOLUTION  OF  THE  STABILITY  PROBLEM 

3 . 1  TheriTiodynaraic  Properties  and  Transport  Coerf Icientg 

Before  numerical  bovmdary  layer  integraftion  can  be  performed 
for  either  the  mean  or  disturbance  equations,  values  for  the 
tharmodyna:nic  and  transport  properties  of  the  liquid  and  their 
required  derivatives  must  be  specified  at  each  step.   Con:>istent 
with  assumptions  made  thus  far  that  the  property  variations  are 
independent  of  pressure  (also  for  the  boundary  layer  approximation, 
the  mean  pressure  is  constant  in  the  boundary  layer  anyway) ,  these 
quantities  can  depend  on  position  only  through  values  of  the  mean 
temperature.   For  the  mean  flow,  such  temperature  dependence 
couples  the  momentum  and  energy  equations.   More  specifically,  for 
the  stability  problem,  accurate  information  is  needed  not  only  for 
the  variable  properties  of  the  mean  flow,  but  for  the  disturbance 
quantities  r,  m,  and  k  as  well  (see  equations  (2.20  and  2.35)). 

Since  the  solution  process  for  the  mean  flow  equations 
requires  chanpes  of  tlic  temperature  at  a  particular  physical 
position  both  during  the  iterative  search  for  unknov.'n  wall  boundary 
conditions  (discussed  in  section  3.2)  and  later  v/ith  the  choice  of 
step  size  for  the  converged  solution,  discrete  input  to  a  numerical 
scheme  of  required  fluid  properties  and  derivatives  is  impossible. 
That  is  to  say,  this  information  is  available  in  f.cncral  on]y  by 


42 


43 


Interpolating  tabulated  experiTrental  data  at  each  step.   It 
foUo'^'s,  than,  that  th.e  temperature  derivatives  inust  be  obtained  b}' 
numerical  dif f erentlati.on  of  the  discrete  set  of  data  points,  and 
so  must  be  interpolated  as  well. 

Water  is  chosen  as  the  representative  liquid  for  numerical 
purposes  to  maximize  applicability  of  the  results  of  this  investi- 
gation (in  \.'hich  thermal  fluctuations  are  considered)  to  situations 
occurring  in  nature,  but  primarily  to  enable  comparison  V7ith  the 
results  of  the  fourth  order  stability  analysis  of  Wazzan,  et.al. 
[4].   In  general  to  obtain  the  "best"  information  available, 
values  for  c  ,  p  ,p  ,  and  k  are  taken  from  recently  published 
measurements  and/or  evaluations  of  preexisting  data.   Temperature 
derivatives  of  these  properties  are  obtained  directly  by 
analytically  differentiating  the  empirical  property-temperature 
relationship  proposed  for  each.   Note  that  this  procedure  is  in 
contradistinction  to  first  estimating  derivatives  directly  from 
the  data  (\:ith  the  method  of  cubic  splines,  for  example)  and  then 
smoothing  any  spurious  fluctuations. 

A  more  extensive  discussion  of  the  sources  and  accuracy  of 
input  property  data,  alternate  m-^ans  considered  for  curve-fitting 
tabulated  data  and  obtaining  high  quality  derivatives,  and  graphical 
presentation  of  Uie  properties  and  temperature  derivatives  for 
water  are  given  in  Appendix  B. 
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3 . 2  Mean  Eqiiations 

The  equations  describing  the.   mean  tx;o-ddmensional,  boundary- 
layer  flow  over  a  heated,  flat  plate  are  obtained  directly  from 
equations  (2.20)  as 


^y 


a    F"/  .-i  rF".o 


and  (3-1) 


,'^  / 


Pr 


z 


with  boundary'  conditions  specified  by  equations  (2.21) 


1 


=  o  .-      f'=  f  =  o  .    1-1  = 


and  (3.2) 

1^  -^  OD  :   F  '  —  1  ,  Hi  —  O 

The  numerical  solution  to  these  equations  introduced  several 
difficulties  characteristic  of  a  boundary  value  problem,  in 
addition  to  those  associated  with  the  nonlinearity  and  coupling 
through  fluid  property  variation.   The  method  of  Nachtsheim  and 
Swipert  [18]  adopted  for  this  analysis  is  found  to  satisfactorily 
deal  v;ith  a]l  difficulties,  especially:  determination  of  v;here 
to  stop  the  integration  for  the  unconverted  solution,  and 
approximation  and  refinement  of  initial  conditions. 

First,  to  obtain  a  unique  solution  of  these  equations. 
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external  conditions  must  be  "satisfied"  for  an  undetermined,  but 
finite  n  =n   _  >  rir.   Asymptotic  converge-nce  (v/here,  F'  -»•  1  and 

H  ->-  0  as  ri  ->  «>  becomes  F'  =  1  and  H  =  0  at  n  =1    )  is  assured 

max 

simply  by  requiring  velocity  and  temperature  gradients  vanish  at 

the  "truncated  infinity"  (F"  =  li'  =  0  at  n=  n    ).   Note  from 

max 

equations  (3.1)  that  these  additional  conditions  cause  all  higher 

derivatives  to  vanish  there  as  V7ell.   However,  the  value  of  n 

max 

(or  equivalently,  that  n  beyond  which  no  appreciable  change  is 

observed  in  the  solution)  still  remains  unknovm.   If  the  a  priori 

specification  of  n  ^  .  is  too  small,  boundary  layer  integration  is 

terminated  too  soon  and  boundary  conditions  cannot  be  satisfied; 

if  T\  is  chosen  too  large,  cormutatlon  time  becomes  excessive, 

max 

Secondly,  numerical  integration  of  the  mean  flow  equations 

requires  specification  of  as  many  conditions  at  the  wall  as  there 

are  conditions  to  be  satisfied  at  the  boundary  layer  edge.   Since 

the  required  additional  wall  conditions  (F"  and  H')  are  initially 

WW  -^ 

unknown,  they  must  be  first  estimated  and  then  iteratively  adjusted 
until  the  asymptotic  conditions  can  be  satisfied,  or 

F'(F"  ,  H'  ,  n    )  =  1 
w    w    max 

H(F",  H';  r,   )  =  0  (3.3) 

w   V7   max 

with  the  supplemental  conditions 

r"(F",  H';  n   )  -  0 
w   v;   max 

H'(F",  H';  n    )  -  0 
w   w   max 

In  other  v;ords,  the  system  of  equations  is  solved  as  an  initial 
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value  problem  with  the  wall  gradients  as  the  specified  initial 
values.   Due  to  the  extrerae  sensitivity  of  equations  (3.1)  to 
values  of  F"  and  H'   even  fair  estimates  of  these  parameters  arc 

W  IT 

likely  to  cause  the  nuraerical  boundary  layer  integration  to  "blow 

up"  prior  to  reaching  the  prescribed  value  of  r\        .   Thus  no  meinis 

would  be  available  for  refining  initial  estimates  of  the  wall 

gradients. 

The  mean  equations  are  thus  solved  in  the  following  fashion. 

VJith  an  estimate  of  initial  conditionn,  F"  and  11"   Runee-Kutta  and 

w      v; 

Adams-Moulton  integration  techniques  are  applied  over  an  abbreviated 
1)   range.   At  the  end  of  this  range,  the  initial  conditions  are 
refined  by  determining  the  least-square  solution  for  discrepancies 
between  the  computed  variables  and  the  proper  asymptotic  values. 
When  successive  iterations  over  this  interval  fail  to  apprecic-bly 
change  the  initial  conditions,  an  error  estimate  E,  the  sum  of  the 
squares  of  the  deviation  of  computed  quantities  from  their 
asymptotic  values,  is  evaluated: 

E  =  (1  -  F')^  +  H^  +  F"^  +  H'^  (3.4) 

If  this  error  is  sufficiently  small,  then  the  solution  is  considered 
to  be  found.   If  such  is  not  the  case,  then  the  ri  range  is  increased 
and  the  above  procedure  repeated.   The  boundary  layer  is  thus 
traversed  in  a  stepv;ise  fashion  until  errors  are  within  acceptable 

limits.   The  last  n  range  then  determines  p    automatically. 

^  max 

In  summary,  this  method  yielded  relatively  rapid,  unitiua 
convergence  to  the  proper  solution  and  seemed  to  be  insensitive  to 
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initial  estimates  of  the  wall  gradients. 

Once  solution  of  the  mean  flow  equations  is  accomplished, 
the  variable  coefficients  for  the  disturbance  equations  can  be 
cctnputed  directly  and  the  latter  solved  in  the  manner  described  in 
the  ne>:t  section.   However,  other  incidental  calculations  can 
easily  be  performed  as  well.   For  example,  to  find  n-''  for  the 
unheated  case  (or  m.ore  specifically,  when  density  is  considered 


constant),  from  equation  (2.18) 


For  the  heated  case,  though,  the  solution  is  not  immediately 
available  from  variables  computed  in  the  boundary  layer  integration: 


f    ^'^         -  F('].nJ  (3.4) 

0  \ 


Instead,  using  stored  values  of  p  (ri)  aiid  the  composite  Simpson's 

Rule  formula  [19,  p. 79],  the  required  integral  is  determined 

num.erically .   Computed  values  of  t\*   are  tabulated  in  Table  E.3  and 

plotted  in  Figure  C.2.   An  analysis  in  section  E.l  indicates  an 

n    =  10  to  be  adequate  for  the  boundary  layer  intejiration  of  the 
max 

mean  equations. 
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Conversion  back    to  j^liysixal   cooidinates   tnust;   be   similarly 
computed.      From  equations    (2.]Ga)    and    (2.17) 


r  - 

7s*  i    -^ 

so   that 

(3.5) 


(3.6) 


a  constant,  once  the  mean  flow  equations  are  solved  and  the 

■*        < 
boundary  layer  edje  is  specified.   A  plot  of  the  X:h  and  —^ 


relationships  for  the  mean  flows  analyzed  is  given  in  Figure  C.l. 

3 . 3  Disturbance  Equations 

The  equations  describing  the  hydrodynamic  spatial  stability 
of  tvjo-dimensionai  disturbances  in  a  two-dimensional,  heated, 
liquid  boundary  layer  are  obtained  from  equations  (2.37-2.40)  by 
setting  6  =  0  and  ;^  =  0  (see  equations  (2.29  and  2.30)),  and 
requiring  that  Re  and  w  =  o<c  be  real  and  c<  be  complex.   Thus, 

Continuity: 
Moraentun: 
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^P 


Energy: 

(u-c)- 


)^M«^<'J  -skp.fe'^^^^^^^'^'t^"^'] 


which  satisfy  a  set   of  hor'oganous  boundary  conditions  obtained  from 
equations  (2.A3-2.45), 


o 


-f  ^  ^  =  r  -^  o 


and 


(3.7) 


1 


oo 


-^  ^,  f  — *-  0 


Since  these  equations  and  boundary  conditions  are  insufficient  to 
establish  a  nontrivial  solution,  the  system  is  cast  as  an 
eigenvalue  problem  (i.e.,  nonzero  solutions  satisfying  boundary 
conditions  are  obtainable  only  for  selected  combinations  of  the 
parameters  a,  uj,  and  Re). 

Since  the  complexity  of  the  equations  precludes  exact 
analytical  solution'^  vjithin  the  boundary  layer,  numerical  techniques 
are  applied.   Of  numerous  possible  methods  previously  utilized  for 
solving  this  type  of  problem,  the  orthonormalized,  differential 
method  used  for  this  investigation  seems  to  be  one  of  the  most 
efficient  and  general  [20].   Schematically,  the  employed  scheme  has 
the  following  steps: 
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(1)  estimate  an  eigenvalue: 

(2)  starting  from  solutions  valid  in  external  region, 
integrate  numerically  to  the  wall  and  tx^y  to  satisfy 
the  specified  boundary  conditions  there  as  veil; 

(3)  refine  the  initial  eij^^envalue  estimate  and  reinte- 
grate as  in  (2) , 

The  final  two  steps  are  repeated  until  "satisfaction"  of  wall 

boundary  conditions  is  attained.  Although  this  procedure  is 

conceptually  quite  straightforward,  its  impleraentation  is  not.   The 

ensuing  discussion  elaborates  on  the  difficulties  associated  with 

each  of  the  steps  enumerated  above. 

As  noted  previously,  there  are  four  parameters  associated 

with  this  eigenvalue  problem,  of  ,  o<  ,  w  ,  and  Re,   In  subsequent 

R    1 

discussion,  "the  eigenvalue"  should  be  interpreted  as  those  tv7o 
real  parameters  which  are  iteratively  adjusted  to  satisfy  the  real 
and  imaginary  parts  of  a  wall  boundary  condition.   The  choice  of 
the  tv7o  fixed  parameters  is  determined  by  information  sought  from 
the  analysis;  for  example,  to  obtain  curves  of  constant  amplifi- 
cation, it  is  found  convenient  to  fix  of  and  Re, and  then  search 
for  the  proper  values  of  c<   and  w  .   Although  inclusion  of  a 
fluctuating  temperature  field  for  this  investigation  is  expected 
to  somewhat  modify  the  existing  results  of  Wazzan,  Okaraura,  and 
Smith  [4]  (in  which  no  such  field  was  considered) ,  the  initial 
eigenvalue  estimate  for  a  specified  T  and  T  is  taken  from  the 
latter.   Once  several  "exact"  eigenvalues  are  determined,  these 
can  be  extrapolated  to  estimate  others  in  the  same  neighborhood. 
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With  all  pertinent  parameters  specified,  the  eigenvalue 
problem  ic  now  converted  to  a  boundary  value  problem.   Thus  it 
remains  to  mathematically  describe  the  asymptotic  condition 
specified  for  the  boundary  layer  edr'.e.   Since  the  mean  flow  co- 
efficients in  the  disturbance  equations  are  all  assuned  to  be 
constant  outside  the  boundary  layer,  an  exact  solution  can  be 


obtained  for  n  >  n  .   Specifically,  equations  (2.42)  reduce  to 


2.-  -^.-^vd 


d; 


0-c)Z, 


K^- 

^5 

3 

Z4 

z;- 

—        L 

2 

Cs' 

l-c)  -f- 


0/ 

Re 


7    - 


to' 


1     Rp  ^5 


(3.8) 


'HAm^-^n. 


e 
— ,  2 


Z,  =   [csRid^c)  +  a-'J  Z,  -isYjf?)^'-^^^3-^  '""^^^ 


Z'  =  [5'+  l3  l?e  Pr  O- C)]  Z 


or  after  some  manipulation 


z,'  -  (5%/)  z,  +-  o^'^'Z,  + 


■1 


.a]o-c)[r^T7;]=o 


dH 


(3.9) 


and 


z;  -  ^  z,  ^-  o 


(3.10) 


v/here 

(3.11) 


S^  --    ^''   +  t-S?2ePr(l-c) 


From  Eqn.  (3.10)  it  is  easily  shown  that  to  satisfy'  the  external 
boundary  conditions  (equation  (3,7)  ),  for  the  j   solution  vector 

Z  J  "  Cj  e  W  "^^  for  n  >  n   and  real  p  >  0 

(3.12) 

Thus  equation  (3.9)  is  an  ordinary,  nonhomogeneous  differential 
equation  with  constant  coefficients  which  when  satisfying  the  same 
conditions  of  boundedness  has  the  solution 

(3.13) 
for  n  >  ri  and  real  a,  y  ^  0« 
From  equation  (3.8)  and  (3.13) 


f--'\^M/e 


for  n  >  n  and  real  a,  y  >  0- 
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The  values  for  Z,  ,  Z  ,  and  Z^_  can  be  exjiresscd  in  terms  of  Z^,  Z^, 


Z     as 

z;^-'  ^ 

^4  0^ 

Re 

T^+^j  v?iH/^ 


(3.15) 


?;i'  =  - « Z3'^' 


As  indicated,  the  three  acceptable,  linearly  independent 
fundamental  solutions  e    '     ,   e     '  '     ,  and   C  ^   (   '  ' 
can  be  combined  to  obtain  a  general  solution  to  the  simplified 
disturbance  equations  outside  the  boundary  layer. 

With  an  eigenvalue  estimate  and  specified  values  for  the 

three  linearly  independent  solution  vectors  Z    at  the  boundary 

(1) 
layer  edge  (obtained  from  eqns.  (3.12-3.15)  by  setting  C,    ~  ^  • -^ 

i,  j  =  1,2,3),  a  fourth  order  Runge-Kutta  method  is  implemented 

to  numerically  integrate  simultaneously  the  disturbance  equation 

solution  vectors  to  the  v?all,  n  =  0.   After  each  integration  step, 

hov/ever,  the  original  linear  independence  is  found  to  be  so  greatly 

diminished  that  the  solution  vectors  become  computationally 

dependent  before  reaching  the  wall.   Thus,  insufficient  information 

is  available  there  to  satisfy  boundary  conditions. 
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In  complex  vector  space,  this  can  be  viewed  as  a  rotation  of  one  or 
more  cf  the  solution  vectors  from  an  orthogonal  to  a  colinear 
configuration,  or  a  reduction  of  tlie  included  angle  from  ninety  to 
•zero  degrees.   The  difficulty  is  explained  as  follows.   Although 
the  exponentially  gro\;ing  solutions  to  equations  (3.16)  can  be 
excluded  from  the  boundary  conditions  in  a  raathematically  exact 
way,  their  influence  cannot  be  eliminated  frora  the  numerical 
integration  inside  the  boundary  layer.  Here,  truncation  or  round- 
off errors  provide  suitable  initial  conditions  for  subsequent 
integration  steps.   As  the  computer  will  allow  the  most  general 
solution,  tlie  rapidly  grov/ing  portion  eventually  dominates  the 
three  solution  vectors  (parasitic  error).   To  avoid  this  difficulty, 
an  orthoromal  basis  is  reestablished  whenever  complete  loss  of 
linear  independence  is  imminent.   Gersting  and  Jankowski  [20,21] 
contrast  this  "near-orthonormalized  integration"  to  the  "fully- 
orthonormalised  integration"  used  by  Wazzan,  et ,  al.  [3],  in  which 
solution  vectors  are  kept  orthonormal  at  each  step. 

Integration  thus  proceeds  in  the  following  manner.   The 
three  basic  solution  vectors  Z    are  integrated  using  a  Runge- 
Kutta  scheme  v/ith  a  fixed  step  size  from  the  boundary  layer  edge 
toward  the  wall  until,  at  mesh  point  n, ,  the  attendant  parasitic 
error  reduces  the  -^ngle  between  any  two  solution  vectors  below  a 
specified  limit,  or 


rr\\n       ^^'~ 


{i''\r^' 


\{t:i'^^ir^\t^^)) 


<   SI         ^J-  1,2.^ 


(3. 10) 

^1 
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where  (Z  "  ,  Z  ^  )  repTesents  th(5  ccnnplcx  irmar  product  of  Lhe 
vectors  1  and  j.   There,  the  Gram-Schmidt  algorithm  for 
orthonorinali^irg  a  set  of  vectors  is  applied  [22],   UTieu 
eigenf unctions  are  to  ba  recovered,  the  ortbonormalizing  matrix, 
available  as  a  nonessential  by-product  of  the  orthonorrnalization 
process,  is  also  coraputed  and  stored  for  later  use.   The  integra- 
tion is  then  continued  from  n-,  using  as  initial  conditions  the 
orthonormalised  basis  vectors,  Z   ,  until  at  iv,esh  point  n^,  the 
angle  between  any  two  solution  vectors  is  again  less  than  fi.   This 
procedure  is  repeated  until  the  v/all  is  reached,  at  vrhich  point 
a  final  orthonorrnalization  is  performed. 

An  evaluation  of  the  proper  step  size,  h.,  and  orthonorrnali- 
zation angle  criterion,  fi,  required  for  this  integration  is  given 
in  sections  E.2  and  E.3  of  Appendix  E  respectively.   Based  on  that 
evaluation,  the  values  of  h,  <  .02  and  ^   =  45''  v.'ere  used  to  generate 

Q 

the  results  presented  in  Chapter  IV.   (Since  the  integration  scheme 
used  requires  computation  of  derivatives  at  half-step  intervals, 
the  mean  equations  thus  must  be  integrated  vjitli  step  size  h  = 
1/2  h^). 

At  the  wall,  the  three  independent  solution  vectors  are 
linearly  combined  in  an  attempt  to  satisfy  the  homogeneous  boundary 
conditions,  equations  (3.7),  for  that  particular  eigenvalue  estimate. 
In  general,  such  satisfaction  is  not  achieved  iiaTiediately  unless  the 
eigenvalue  is  "exact,"  and  so  some  iterative  correction  scheme  must 
be  implemented.   Thus,  for  the  vector  components 
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Z.(o)-^CjZ/(o) 

J"  I 

it  is  required  that  (see  also  [23]) 


(3.17) 


Z  ^o)  4-  B,  2ao)  =0 


and  (3.18) 

where  specification  of  the  arbitrarj-  constant  B  ~  u   (0)  is  found 
to  give  good  eigenvalue  convergence  using  the  sesrch  procedures 
discussed  below.   Following  Wazzan,  et.  al.  [23],  the  quantity 
actually  used  to  clearly  distinguish  an  eigensolution  from  all 
others  is  the  complex,  surface-normal,  boundary  layer  adrrdttauce 


^^(o) 


for,  the  smaller  the  "zero"  components  of  the  eigensolution 

(Z  (o) ,  Z_(o),  and  Z„(o))  are  than  the  nonzero  ones  (Z^  (o) ,Z^ (o) ,  and 

Z,  (o)),  the  more  accurate  will  be  the  eif.envalue  [24,  p.  280]. 

D 

Requiring  both  the  real  and  imaginary  parts  of  Y  vanish  establishes 
two  real  equations  in-  the  two  adjustable  parameters  (or  the  complex 
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"eigenvalue") . 

For  iterative  refinement  of  the  eigenvalues  (i.e. ,  finding 
the  2ero(s)  of  Y  )  ,  tv.'o  schemes  ware  employed:  the  plane  fit  of 
Wazzan,  et.al  [23],  and  the  quadratic  complex  curve  fit  of 
Muller  [25],   In  the  former,  a  plane  of  the  form 

"eigenvalue"  =  a  +  aA'  _  +  a„Y  ^ 
'^  o    1  oR    2  oi 

is  fitted  to  three  successive  eigenvalue  estimates  and  their 

corresponding  computed  admittances.   Tlie  corrected  eigenvalue 

estimate  for  the  next  iteration  is  then  simply  a  .  \Jhen   Y  is 

*^  -^   o         o 

sufficiently  snail,  Muller 's  method  is  used  to  expedite  eigen- 
value convergence.   This  scheme  fits  a  quadratic  curve  in  complex 
space,  also  through  the  lust  three  eigenvalue  estimates  and  com- 
puted admittances,  of  the  form 

2 

Y^  =  b  +  b^  (eigenvalue)  +  b_  (eigenvalue) 

Setting  Y  =0,  this  equation  can  be  inverted  using  the  standard 

quadratic  formula  to  give  an  eigenvalue  estimate  for  the  next 

iteration.   Note  that  both  techniques  have  the  advantage  that  no 

derivative  of  Y  and  only  one  evaluation  of  Y  is  required  per 
o  o 

iteration  (in  contrast,  say,  to  Newton's  method).   A  complete 
discussion  of  the  more  subtle  applications  of  these  tv;o  procedures 
for  an  efficient  eigenvalue  search  is  relegated  to  section  E.4. 

VJhen  the  "eigenvalue"  is  determined  to  v.'ithin  desired 
accuracy,  the  eigenfunctions  can  be  recovered  using  the  stored 
orthonormalizing  matrix  and  solution  vectors  [22].   For  standardi- 
zation, the  magnitude  of  these  eigenfunctions  arc  adjusted  so  as  to 
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make  the  amplitude  of  the  pressure  fluctuation  equal  to  unity  at 
the  v.'all.   If  desired,  energy  production,  dissipation  and  transfer 
terms  can  then  be  evaluated  directly  (see  section  2.4). 


ClL/vPTER  IV 
RESULTS 

The  hydro dynamic  stability  of  the  boundary  layer  formed  on 
a  heated,  constant-temperature  flat  plr.te  (M  =  0)  in  tX'Zo-dimensional 
water  flow  and  with  tuo-dimensional  disturbances  {^  -  0   and  6=0) 
is  examined  and  the  results  presented  in  this  section.   Specific 
comparison  is  sought  with  the  results  of  VJazzan,  et.al  [4]  vjho 
analyzed  this  same  configuration,  but  who  assumed  only  isothermal 
disturbances  of  the  mean  variables. 

Rather  than  regenerating  in  this  investigation  the  stability 
characteristics  for  the  entire  range  of  wall  temperatures  considered 
by  Wazzan,  et.al,,  selective  checks  are  made  to  compare  the  results 
of  the  tv.'o.   Accordingly,  since  those  authors  found  that  increased 
surface  heating  at  first  stabilized  and  then  destabilized  the 
boundary  layer  flow  (i.e.,  peaking  of  Re  .     .    -  that  value  of 
Re  below  which  disturbances  are  damped  for  all  wave  numbers  and 
frequencies),  V7ith  the  extremums  occurring  at  about  T_  =  140°F 
(T  =  60°F)  ,  the  numcric.-^l  calculations  herein  will  be  limited 
primarily  to  wall  temperatures  v;ithin  the  nomial  liquid  range  of 
water  and  including  that  temperature  range  where  the  peaking  occurs: 
specifically,  T  =90°,  ]50%  and  200°F  with  T^  =  60°F.   V;hen 
detailed  characteristic  behavior  is  to  be  illustrated,  the  case  of 
T  =  90°F  will  be  used,  as   this  is  conceivably  within  the  range  of 
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temperatures  for  which  e>qperiiiiental  comparison  can  be  made. 

An  accurate  solution  to  the  mean  fiov;  equations  is  essential 
for  solving  the  disturbance  equations.   The  procedure  and  results 
for  the  mean  flow  solution  are  described  in  Section  3.2  and 
Appendix  C.   In  the  latter  both  detailed  boundary  layer  distri- 
butions of  all  the  mean  properties  and  their  derivatives  required 
for  the  coefficients  of  t'ne  disturbance  equations ,  and  variations 
of  surface  velocity  and  temperature  gradients  with  the  wall 
temperature  are  graphically  presented. 

The  ir.ean  velocity  and  temperature  boundary  layer  profiles  to 
be  examined  in  this  section  are  repeated  in  Figure  4.1.   Note  that 
the  temperature  boundary  layer  is  much  thinner  than  that  of  the 
velocity  boundary  layer  so  thermal  influences  are  restricted  to 
regions  close  to  the  v;all.   This  observation  is  related  to  the  later 
understanding  of  the  relative  importance  of  thermal  disturbances  in 
boundary  layer  stability  calculations. 

Keeping  in  mind  the  role  that  a  linear  stability  analysis 
should  play  in  only  suggesting  those  conditions  which  might  enhance 
or  delay  the  onset  of  laminar-turbulent  transition  (there  is  no 
definitive  link  between  the  point  of  instability  predicted  by  such 
an  analysis  and  the  point  of  transition) ,  this  investigation  is 
restricted  to  regions  which  could  conceivably  still  be  laminar. 
Thus,  although  determination  of  stability  characteristics  for 
excessively  large  Re  might  bo  an  interesting  numerical  or  mathemati- 
cal exercise,  tlie  results  are  not  of  interest  if  transition  has 
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occurrea. 

Jluch  of  this  investigation  was  conducted  neglecting  the 
effects  of  a  fluctuating  density  field  (r  =  0);  that  is,  only 
fluctuations  of  viscosity  in(y)  and  thermal  conductivity  K(y)  were 
included  in  the-  disturbance  nion;entura  and  energy  equations,  respec- 
tively.  Such  an  assuraption  seemed  reasonable  due  to  the  relative 
constancy  of  the  density  of  water  over  the  normal  liquid  range. 
Subsequent  calculations  with  variable  density  nevertheless  were 
made  and  their  results  are  included. 

Unless  other^;ise  indicated,  the  propertji'-teTTiperature  varia- 
tion for  v.'ater  is  that  specified  in  Appendix  B.2.   Although  this 
variation  differed  slightly  from  that  of  Kaups  and  Smith  [7]  (used 
by  Wazzan,  et.al.  in  their  analysis)  as  shovzn  in  Figure  B.7,  there 
is  a  marked  effect  on  some  of  the  calculated  stability  character- 
istics, as  will  be  shown  presently. 

4  . 1  Stability  Characteristics 

As  suggested  in  previous  discussion,  for  a  specified  set 

of  mean  flow  parameters  (T  ,  T  ,  and  M) ,  the  eigensolutions  to  the 

homogeneous  disturbance  equations  and  their  boundary  conditions 

describe  a  surface  in  the  four-diirensional  parameter  space 

(a„ ,   a^  ,  u   ,  and  Re).   It  is  of  interest  to  locate  minina  and 
R    I    " 

maxima  on  this  surface  and  to  trace  the  path  followed  by  an  actual 
disturbance  as  it  propagates  through  the  boundary  layer  (e.g.,  to 
ascertain  v;hcther  an  arr.plif ication  region  (a  <  0)  is  traversed, 


and  if  so,  the  rate  of  amplification  expected).   The  stability 
characteristics  associated  with  neutral  and  amplified  disturbances 
are  discussed  below. 


4.1.1  Neutral  Disturbances  (a  =  0) 


Perhaps  the  most  obvious  quantitative  guide  to  the  effect  of 

heating  on  boundary  layer  stabilization  is  the  amount  vith  which 

neutral  curves  are  shifted  with  changes  in  the  wall-to-free  stream 

temperature  difference.   The  movement  of  these  stability  loops  is 

characterized  by  changes  in  the  minim.um  critical  Reynolds  number. 

Re  .     ,  ,  the  maxiiaum  wave  number  a^    ,  and  the  maximum  fre- 
min.crlt.'  R 

max 

quency,  to     (above  which  all  disturbances  are  damped) .   As  may 

max 
be  seen  from  Figures  4.2  and  4.3,  with  only  moderate  licating,  both 

the  fourth  order  sj'stem  of  equations  of  Wap^zan,  et.al.  [4]  and  the 
sixth  order  system  of  tliis  investigation  (vjithout  density  fluctua- 
tions, however)  predict  a  significant  delay  in  the  onset  of 

boundary  layer  instability  (evidenced  by  the  increased  Re  .     ,   ) 
■'        ■'  ■'  mm.crit. 

and  a  reduced  range  of  frequencies  and  \-iave   numbers  which  could 

become  unstable.   As  the  vjall-to-f ree  stream  temperature  difference 

is  further  increased,  however,  both  the  fourth  and  sixth  order 

systems  consistently  indicate  a  peaking  and  subsequent  decrease  of 

Re  .     .   .   Based  on  these  calculations  for  water  at  the  indicated 
mm.  crit . 

temperature  levels,  the  following  conclusions  may  be  dravm. 

(1)  It  would  be  impossible  to  achieve  cor-.plete  boundary 
layer  stabilization  by  indefinitely  increasing  the  heat  rate; 
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O   V/azzan,  Okaraura,  and  Smith  [A] 
□   Present  investigation  (i"~0) 

All  points  are  computed  using  the 
property-temperature  variations  of 
Kaups  and  Smith  [7j. 


-p-1 — I — r— r"r— r 

T-r-T-T 

300 

T^    Cf) 

Fipura   4.3 


Con-ir>arison  of  Re    .  ,         computed  with    and  v/tthout 

min.crit. 


thermal  disturbances   as   a   function   of  T      (fixed  T 

w  " 


50"f) 
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(2)  ifi  lie  the  calculations  \;!rJ.ch  inv-,i  i^;de  thermal  disturbances 

do  predict  slic^htlv  larger  values  of  Re  .     .   ,  the  qur]  itati  ve 
^  mm.  crit . 

variation  of  this  quantity  v/ith  T   (for  fixed  T  =  60"?)  is  not 

altered. 

The  values  of  Re  .     .'       are  easily  obtained  for  specified 
Kin.crit. 

mean  flow  parameters  by  fixing  a^  =  0  and  a  »  than  adjusting  esti- 
mates  of  Ra  and  (o  until  boundary  conditions  can  be  satisfied .   An 

CO  •' 

example  of  this  procedure  is  shown  in  Figure  4,4  and  4.5  for  several 
different  ca;;e.-:. 

For  no  other  calculation  is  the  influence  of  the  chosen  pro- 
perty-temperature relationship  or  assumed  property  fluctuations  more 
pronounced  than  in  this  one.   For  example,  from  these  two  figures 
it  may  be  seen  that: 

(1)  use  of  the  property  variation  of  Kaups  and  Smith  {7] 

effects  a  shift  of  Re  .     .    first  above  (T  -  90°F) 

mm.crit.  v; 

and  then  below  (T  ^  --  200° F)  that  computed  vrith  the 
variation  specified  in  Appendix  B.2; 

(2)  for  either  property  variation,  calculotions  including 

thermal  disturbances  (but  with  r  =  0)  predict  a 

Re  .     .    roughly  10%  greater  than  that  computed 
mm.  crit . 

using  the  fourth  order  equation  of  VJazzan,  et.al. 
(equation  (1.2)) ; 

(3)  inclusion  of  density  fluctuations  in  the  calculations 
as  vi;ll,  offsets  the  gains  in  stability  predicted  above 

in  (2)  as  T  is  increased. 

w 

An  exp.1anation  for  the  relatively  pronounced  infl>;ence  of 

density  fluctuations  on  Re  .     .    calculations  is  found  in  the 
■'  mm.crit. 

distribution  of  disturbance  properties  in  the  stability  equations. 
Although  the  density  fluctuations  may  be  comparatively  small  (lay 
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as  much  as  several  orders  of  magnitude  less  than  the  viscosity 
fluctuations  -  see  Figures  4.9  and  A. 10,  for  example),  they  do 
appear  in  both  the  inviscid  and  viscous  parts  of  the  stability 
equations,  whereas  viscosity  fluctuations  are  restricted  only  to 
the  latter.   As  temperatures  are  increased,  the  disparity  between 
the  two  disturbance  quantities  diminishes;  eventually,  disturbance 
density  contributions  dominate  those  from  all  other  property 
fluctuations  in  the  calculations. 

Hov7ever,  as  shown  in  Figures  4.4  and  4.5,  while  the  fluctu- 
ating viscosity  terms  tend  to  enhance  stability,  the  density  terms 
degrade  it.   It  is  important  to  note  that  over  the  norm.al  liquid 
range  of  water,  results  from  the  sixth  order  system  of  equations 
(even  including  density  disturbances)  predict  a  slightly  more  stable 
boundary  layer  than  does  the  fourth  order  system  for  which  r  (y)  =  0. 

4,1.2  Amplified  Disturbances  (a^  <  0) 

The  tendency  to  transition  is  also  a  function  of  amplifica- 
tion rates  and  so  attention  will  now  be  given  to  the  amplified 
portion  of  the  stability  loop  (a  <  0). 

A  typical  set  of  stability  curves  illustrating  the  topologi- 
cal nature  of  the  eigenvalue  variation  for  T  ~-   90°F,  T  -  60°F,  and 
M  =  0  (r  =  0)  is  shown  in  Figure  4.6.   These  curves  are  similar  in 
general  shape  to  those  for  the  unheated  Blasius  boundary  layer 
although  dramatically  different  in  magnitude  (see  reference  [ Z6 ,  p. 90] 
for  example) . 
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A  much  more  sicnifj.cant  indication  of  boundary  layer  stabi- 
lization by  heatinn  than  Re  .     ,    is  eiven  by  the  inaKimum 
■^        -        inm.crit.  ^ 

amplification  rate,  -^cx  /Re)   ,  established  for  specified  mean  flov? 

conditions.   The  paraiaetric  location  at  which  this  maximum  occurs 

is  determined  by  a  more  involved  procedure  tb.an  that  for  finding 

Re  .     .   ,  as  illustrated  in  Figure  4.7.  More  specifically,  the 
rain.crit. '  *-  r  }  ■, 

amplification  rates  for  specified  Re  are  computed  and  individual 
local  m.axima  determined.   Follov/xng  the  locus  of  these  individual 
maxima,  the  largest  value  attained  is  designated  the  maximum  ampli- 
fication rate  and  its  associated  wave  number  (and  frequency)  are 
obtained  graphically.   IsTiile  the  flatness  of  this  locus  curve  in  the 

neighborhood  of-fa-r/Re)     allov/s  relatively  accurate  evaluation 
\  \       ■'  max 

of  that  quantity,  it  does,  however,  preclude  similar  accuracy  for 
Op  (and  u^)  there.   Physically  this  implies  that  there  is  effectively 
a  band  of  v/ave  numl^ers  (and  frequencies)  V7hich  will  receive  nearly 
maximum  amplification.   To  complete  the  eigenvalue  set,  Re  can  be 
determined  from.  Figure  A. 6  by  again  following  the  locus  of  indivi- 
dual maxima  until  the  required  value  of  a  is  reached. 

Due  to  the  numerical  effort  required  to  compute  the  maximum 
amplification  rate,  only  one  set  of  calculations  v/as  m.ade  for 

T  =  90°F  and  T  =  200°F  (T   =  60°F) ,  using  the  presumably  more 
v;  V7  <^ 

accurate  property-temperature  relationships  of  Appendix  B .  2  (assum.ing 
r  =  0).   Values  obtained  for  -(a^/^'c)  ^.^^  x  lO'  of  26.6  for  the  former 
and  19.5  for  the  latter  can  be  compared  with  the  corresponding 
values  tabulated  by  Uazxan,  et.al.  [4]  of  28  and  20,  and  the  value 
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for  the  unheated  (Blnsiuc)  boundarj'  of  7^iO. 

Thus,  not  on]y  is  the  onset  of  disturbance  instability  deJayed 
and  the  range  of  disturbance  wave  numbers  and  frequencies  receiving 
amplification  reduced  as  pointed  out  in  tb.e  preceeding  section,  but 
the  ariplif ication  rates  to  which  disturbances  v;ill  be  subjected  are 
also  draruatically  reduced.   These  features  all  point  to  a  Piore 
stable  boundary'  layer  flow. 

4 , 2  Disturbance  Amp litude  and  Energy  Distributions 

Finding  no  major  modifications  to  the  stability  characteris- 
tics of  a  heated,  V7ater  boundary  layer  by  including  in  the  calcula- 
tions the  effects  of  a  disturbance  temperature  field,  an  explanation 
is  sought  from  the  boundary  layer  distribution  of  disturbance 
quantities  and  the  energies  produced  or  dissipated  by  them. 

Due  to  the  hoinogeneit)'  of  the  equations  and  boundaTry  con- 
ditions, solutions  to  the  disturbance  equations  are  independent  of 
an  absolute  amplitude  (ss  section  2.3.1).   Thus,  as  a  comparison 
is  sought  betA<7een  disturbance  levels  for  both  the  same  and  different 
mean  thermal  boundary  conditions,  it  is  necessary  to  resort  to  some 
physical  reasoning  to  determine  a  normalizing  quantity  that  v;ould 
remain  virtually  constant  in  the  comparison.   Since  there  is  no 
strearcwise  pressure  gradient  (i.e.,  M  =  0  for  all  Re),  tlic  mean 
pressure  v/lthin  the  boundary  layer  is  assumed  constant,  and  the 
disturbance  pressure  is  nonzero  and  slowly  varying  throughout 
(see  Figures  4.8  through  4.10,  for  example),  the  amplitude  of  tlie 
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Fijiuro   4.9     Neutral   cii.sturb;inct;   amplitude  profile,';   for   the  lieatcd 
boundary   layer:   T     °   9o"f   and  T^  =   60°F. 
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disturbance  pressure  at  the  wa.l]  ,  [rr  (0)  [ ,  vas  chosen  to  normalize 
all  disturbance  quantities.   Thus  all  comparisons  made  betxireen 
different  eigensolutioas  in  the  remainder  of  this  section  are  pre- 
dicated on  this  assumption. 

To  further  investigate  the  importance,  of  density  fluctua- 
tions, they  are  included  in  all  calculations  made  in  this  section. 

Typical  amplitudes  of  disturbance  quantities  found  in  this 
analysis  are  presented  in  Figures  4.8  through  4.10  for  wall  tempera- 
tures T^  =  60°,  90%  and  200°F  (T,^  =  50"?),  respectively.   For 
standardization,  the  eigenvalues  associated  with  each  set  of  eigen- 
functions  are  all  located  on  the  upper  branch  of  their  respective 
neutral  stability  curves  at  frequencies  near  that  receiving  maximum 
amplification. 

After  studying  these  figures,  the  follovring  observations  are 
made. 

(1)  As  anticipated  from  the  relatively  large  Pr  for  water 
(i.e.,  such  that  the  thermal  boundary  layer  thickness  is  much  less 
than  the  velocity  one),  in  the  temperature  range  considered,  thermal 
effects  are  restricted  to  a  region  very  close  to  the  \;all.   Since 
theriual  fluctuations  must  vanish  at  the  wall,  even  within  this 
region,  they  attain  fairly  large  ara-plitudes  only  within  a  far  narrov/- 
er  range. 

(2)  By  far  the  most  dominant  property  disturbance  amplitude 
is  that  of  viscosity.  However,  for  relatively  small  surface  heating 
(T  =  90*'F),  even  the  thermal  conductivity  disturbance  amplitudes  is 
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commensurate  with  that  of  the  normal  velocity  in  some  regions. 

For  large  heating  (T  =  200°F)  ,  \<ihen   density  fluctuations  are 

w 

significant,  they  are  seen  to  be  larger  iu  amplitude  than  both  the 
thermal  conductivity  and  normal  velocity  disturbances. 

(3)  With  increased  heatiuf;,  the  ratio  of  f    /  I  p  6 1 

'  '  max   '  "^  ^  '  max 

increases  significantly.   Further,  assuming  that  the  disturbance 

wall  pressure  amplitude  is  relatively  insensitive  to  changes  in 

T  within  the  range  60"  iS  T  ^  200°F,  If  I    increases  monotonically 
w  '^  w        '  '  'max  -^ 

with  T  while  I  p  cj     exhibits  a  m.aximum. 

(4)  Just  as  the  stability  characteristics  are  only  slightly 
modified  by  the  inclusion  of  thcmal  disturbances  in  the  analysis, 
so  too  are  t;he  disturbance  velocity  and  pressure  amplitude  pro- 
files.  Even  the  disparity'  shovrn  (except  in  the  region  where  thermal 
disturbances  obviously  doirdnatc)  can  be  attributed  in  part  to  the 
difference  between  the  neutral  eigenvalues  required  to  satisfy  the 
fourth  and  sixtVi  order  systems. 

Once  the  eigenvalues  and  eigenf unctions  are  obtaJned  by  the 
procedure  outlines  in  section  2.3.1,  it  is  a  simple  matter  to  com- 
pute the  energy  balance  developed  in  section  2. A.   Hov7ever,  for  ease 
of  analysis,  only  those  terms  contributing  a  net  increase  or  de- 
crease of  energy  of  a  neutral  disturbance  are  evaluated  (i.e.,  all 
term>s  of  equation  (?.  .50)).   The  results  of  sucli  calculations  are 
presented  in  Figures  4.11  through  4.16.   For  cacli  of  the  cases 
T  =60%  90°,  200°F  (T  ^  60" F),  the  eigenvalues  selected  from  the 
upper  and  lower  branches  of  the  neutral  stability  curves  (see 
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Figure  4.11   Enerry  discributions  of  a  neiitml  disturbance  In  the 
unheated  (Blasin-o)  boundary  layer  (lower  hrr.nch) . 
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Fij'iirp  4.12   Energy  di.TtrJ.butlons  of  a  neutral  disturbance  fa  the 
unhcatcd  (Blaslus)  boundary  layer  (upper  branch). 
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FJj'ure  4.13   Energy  distributions  of  a  neutral  dit; turbancc  in  the 
heated  boundary  layer:   'J'  =  90  F  and  T  f   00  F  (Jower  branch). 
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Ip.ure   A.  14      Energy   cii.stribuLton  of   a   ncucra].   disturbance    in   the 

heated  boundary    layer:    T     =   90°F   and  T     =   60°F    (upper  branch). 
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Figure  4.15   F.nerj'.y  dltitr ibutions  of  a  neutral  dinturbanco  in  the 


hcatc^i  boundary  Layer:   T  ■■=   200  V   and  'i. 
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CO  F  (lovor  branch). 
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Figure   A. 16        Energy   d) -jtrJlnitlons   of   a   neutral    diyturbance    In    the 


heated   boundary    layer:      T     =■    200"f   and   T      »   60°F    (upper   brancli) . 
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Figure  A.l)  again  are  for  frequencies  for  which  disturbances  are 
subjected  to  approximately  maximuin  arnplifi cation  v;hen.  traversing 
the  loop . 

From  these  six  figures,  the  following  observations  are  made. 
The  "terms"  referred  to  are  those  designated  in  equation  (2.50). 

(1)  The  dominant  terras  in  the  neutral  disturbance  energy 
balance  are  the  traditional  ones:   The  Reynolds 
stress  production  (term  (2))  in  the  region  between 
the  vjall  and  the  neighborhood  of  u  =  C  ,  balanced  by 
the  dissipation  (term  (6))  which  is  extreraely  large 
near  the  v;all. 

(2)  Dissipation  due  to  a  nonunifoi'n  mean  disturbance 
density  (terra  (7))  is  always  negligible.  Although 
not  nearly  as  small,  energy  production  due  to  a  mean 
viscosit}'  gradient  (terra  (5))  is  also  insignificant. 
Neither  are  large  enough  to  appear  in  Figures  4.11  - 
4.16. 

(3)  The  disturbance  shear  stress  (terra  (3))^  although 

dissipative  in  some  regions,  represents  primarily 

energy  production  for  the  disturbance  field.   In 

all  the  cases  examined,  it  produced  roughly  10%  of 

the  energy,  v;ith  the  Rej'nolds  stress  production 

(term  (3))  supplying  the  rest.   Recall  that  this 

is  approximately  the  same  amount  that  Re  . 

mj.n.crxt. 

was  increased  for  sixth  order  system  calculations 

(r  =  0)  beyond  that  calculated  v;ith  the  fourth 

order  system  (see  Figures  4. A  and  4.5). 

(4)  Pressure  production  (term  (4)),  also  resulting 
from  density  nonunjLf ormities ,  primarily  gives_rise 
to  production  for  u  >  c  and  dissipation  for  u  <  c  . 

Since  its  integrated  value  over  the  boundary  layer  is 
nearly  zero,  this  term  is  not  thought  to  be  of 
significance  in  the  stability  of  this  flow. 

To  check  the  accuracy  of  the  computed  energy  production  and 
dissipation  terms,  each  was  first  numerically  integrated  from  the 
wall  to  the  boundary  layer  edge  and  then  the  resultant  summed.   With- 
in the  accuracy  sought,  this  sum  was  zero  as  required  by 
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equation  (?.50).   Furthermore,  it  vos    found  that  the  percent  of 
positive  or  negative  energiy  contribution  from  the  production  terms 
(such  as  the  Rtynolds  stress)  reniained  virtuallj"  unchanged  with 
changes  in  the  v;all  temperature  level. 

It  is  especially  ir.'portant  to  note  the  extent  of  the  region 
in  v/hich  significant  energy  contributions  take  place.   For  the 
heated  boundary  layer,  in  contrast  to  the  unheated  one,  all 
energetics  are  restricted  to  a  region  very  close  to  the  wall.   Out- 
side of  that  region,  energy  terms  tend  to  vanish,  signifying  that 
there  is  very  little  correlation  betv^een  disturbance  quantities  over 
most  of  the  boundary  layer. 


CHAPTER  V 
Sm-QIARY  OF  RESULTS  AND  CONCLUSIONS 


This  investigation  has  Gxarriined  the  sixth  order  system  of 
equations  describing  the  hydrodynamic  stability  of  infinitesimal 
disturbances  in  a  heated,  flat  plate,  water  boundary  layer  flow. 
Throughout,  comparison  has  been  made  with  the  fourth  order 
analysis  of  this  same  problem  as  previously  done  by  Waz^an, 
Okamura,  and  Smith  [4] . 

The  time-averaged  (mean)  equations  of  laminar  flow  are 
simplified  by  making  the  usual  boundary  layer  assumptions,  and  by 
neglecting  buoyancy  forces,  viscous  dissipation  and  expansion 
energy.   They  are  then  transformed  from  partial  to  ordinary 
differential  equations  using  Howarth-Dorodnitsyn  type  similarity 
variables  and  their  resultant  solutions  used  to  provide  the  velocity, 
temperature,  and  other  profiles  required  for  ttie  disturbance  equa- 
tions.  The  latter  are  found  as  the  linearii'.ed  difference  between 
the  instantaneous  and  mean  equation,  after  neglecting  here  as  well 
buoyancy,  dissipation,  and  expansion  energy  terms,  assuming  the  flow 
to  be  locally  parallel,  and  finally  using  a  normal  modes  analysis 
to  obtain  ordinary,  linear  differential  equations.   The  tvjo  sets  of 
equations,  mean  and  disturbance,  thus  derived  differ  from  those  of 
Wazzan,  et.al.  by  the  inclusion  of  a  disturbance  energy  equation 
coup]ed  to  the  other  disturbance  equations.   This  permits  density. 
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viscosity,  and  thermal  conductivity  fluctuation  effects  to  be 
considered  in  the  stability  problem. 

Eoth  the  fourth  order  system  of  Wazzan,  et.al.  and  the  sixth 
order  system  of  the  present  investigation  predict  significant 
boundary'  layer  stabilization  with  moderate  heating,  but  display  a 
maximum  and  subsequeiit  stability  decrease  as  the  plate  surface-to- 
free  stream  ternperature  difference  is  further  increased.   The 
computed  stability  characteristics  (minimum  critical  Reynolds  Number, 
maximum  amplification  rate,  maximum  wave  number  and  frequency  sub- 
ject to  amplification,  etc^  obtained  from  each  system  differ  only 
slightly,  however.   In  particular,  the  peaking  of  each  characteristic 
with  V7all  temperature  (at  about  T^  =  lAO^F  with  fixed  T^  =  60"F)  is 
neither  eliminated  nor  significantly  altered  by  the  higher  order 
system.   Thus,  complete  stabilization  apparently  cannot  be  achieved 
by  indefinitely  increasing  boundary  layer  heating. 

The  insensitivity  of  these  results  to  the  additional  con- 
sideration of  thermal  disturbances  is  attributed  to  the  extremely 
limited  region  in  which  such  disturbances  are  important.   Du2  to  the 
relatively  large  Prandtl  Number  for  v/ater  and  homogeneous  dis- 
turbance boundary  conditions,  the  boundair)'  layer  distributions  of 
thermal  disturbance  amplitudes  "spike"  very  close  to  the  wall. 

In  all  cases  considered,  the  sixth  order  stain  !ity  system  is 
found  to  predict  a  slightly  more  stable  flovj  than  the  fourth  order 
system,  as  evidenced  by  larger  values  of  the  minimrm  critical 
Reynolds  Number  and  reduced  values  of  the  maximum  amplification  rate. 
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However,  it  is  postulated  that  for  wall  temperatures  beyond  the 
normal  liquid  range  of  water  (with  T  =  60°F) ^  the  presence  of 
density  fluctuations  v;ill  eventually  negate  the  increa.s£d  stability 
predicted  with  only  viscosity  fluctuations  assimed,  and  actually 
degrade  the  degree  of  stability  predicted  by  th.e  fourth  order 
system.   That  is,  when  thermal  disturbances  are  assumed,  it  is 
necessary  to  include  all  property  disturbances  to  obtain  an  accurate 
assessment  of  the  boundary  layer  stability. 

The  admission  of  thermal  fluctuations  also  has  revealed  the 
existence  of  several  energy  production  and  dissipation  terms  not 
present  in  an  energy  balance  equation  for  the  fourth  order  system. 
Of  particular  note  is  the  disturbance  shear  stress  production, 


—  i?   ^ +  —v—    Jd-kL  ,    a   source  v^hich  accounts  for  about  10%  of 

the  energy'  extracted  from  the  mean  and  supplied  to  the  disturbance 
flow.  In  contrast,  when  onl}'  isothermal  disturbances  are  assumed, 
energy  is  supplied  completely  by  the  Reynolds  stress  production. 

In  surmiary,  the  fourth  order  system  of  equations  used  by 
Wazzan,  et.al.  [A]    scens  to  adequately  describe  the  stability  of  a 
heated,  v;ater  boundary  layer  flovj  in  the  temperature  ranges 
considered.   This  may  not  be  true,  however,  for  other  liquids, 
particulr.rly  those  vnth  Pr  <  1,  v;here  it  is  conjectured  that  calcula- 
tions will  be  much  more  sensitive  to  the  influence  of  thermal  dis- 
turbances. 
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APPENDIX  A 
Detailed  Derivation  of  Mean  and  Disturbanc e  Flow  Equati o n s 

In  accordance  with  the  discussion  in  section  2.1,  it  can 
be  assumed  that  the  flo^;  field  is  describable  by  the  set  of 
conservation  equations  (2.1-2.3).   Mean  and  fluctuating  behavior 
are  deduced  from  these  equations  in  the  following  fashion. 

The  instantaneous  value  of  each  quantity  can  be  decomposed 
into  a  mean  and  fluctuating  component. 

Q,  =  Q.(Y^)  ^  4U  ,t)  (2.G) 


where  Q.  is  the  time-averaged  mean 


'*-''^  -   -     -  (A.l) 


J  I  ».  CO     I    J  ■' 


and  vjhere  the  time-averaged  value  of  any  fluctuating  component 
is  defined  to  be  xero. 


Q   =Tl»Tr^«'T,,r)»fr=o 


For  a  time  average  to  be  nieaningfuljmean  quantities  must  be 
time-independent  as  already  denoted  by  the  argument  of  0.  in 
equation  (A.l).   The  use  of  time-   rather  than  space-averaged  means 


94 


95 


corresponds  to  observed  boundary  layer  behavior  v:here  disturbances 
grow  or  decay  spatially  rather  than  temporally. 

Thus,  substituting  equations  (A-1)  into  equations  (2.1-2.3) 
and  perfomlnf  a  tine  average,  the  mean  equations  become 
(dj  iTiensionally ) 


Continuity: 


(  0U.   +  o  U  ^  j  =  O  (A. 2) 


Hotnentuvi: 


=  -^- 13,  - 1  ^P+  ^£,_,.  +-  >,?:_) 4-zi  ^;i6„- v^s 


J 
Energy: 


"'^'f  "jsi    =(P"e^"jf^^^^    <? 


'J ij    -^ 


Eventually  (see  section  2.3),  to  reduce  the  disturbance 
differential  equations  from  partial  to  ordinary  ones,  disturbances 
of  tbe  form  3,  =  ?(/;  e^p  [f  i'=^ft-=^iX''xcose -t  Zii'^  fi)-tiOt]  will  be  considered. 
Clearly  the  tinit^  average  of  this  disturbance  is  zero,  whereas  a 
space  averap.e  is  not  (except  for  a  neutral  disturbance  where  o<  =  0). 
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To  avoid  the  closure  prcblem  characteristic  of  turbulent 
flo\.'s,  it  is  desirable  that  this  set  of  equations  be  solved 
independently  of  correlated  disturbance  quantities.   Thus,  a  set  of 
nev7  nondimsnsional  independent  variables  is  defined 


h  '  ('-  t')X 


+ 


-J. 


w-ri-r^ 


v;hare  L.  and  1.  are  lenpth  scales  in  the  ith  direction  character- 
istic  of  mean  and  fluctuating  flow  behavior,  respectively.   For 
example,  at  a  given  position,  L  might  represent  a  length  character- 
istic of  the  boundary  layer  development,  specifically  the  distance 
from  the  body's  leading  edge,  while  1_   v7ould  correspond  to  the 

disturbance  wave  length  at  that  point.   When  characteristic  length 

1. 
scales  are  comparable  so  that  — ' —  =  0(1),  then  equation  (A. 5) 

reduces  to  E,.    -  x..      Here  x,  most  be  considered  a  measure  of  both 
3     3         3 

large  and  small  scale  variations  in  the  j th  direction,  since  both 

scales  are  iu  fact  identical,  or  formally  £.  =  X.  =  x..   Such  a 

3     3     3 

situation  typically  arises  iii  this  analysis  where  transverse 
boundary  layer  scales  L  and  1  are  equivalent  (exception  to  be 
noted  later)  and  equal  to  the  boundary  layer  or  displacement 
thickness,  for  example,  so  that  y  =  Y. 

If  it  can  be  argued  that  variations  in  the  mean  motion  be 
Independent  of:  small  scale  disturbances  in  each  direction  (i.e., 
mean  equations  can  be  solved  independently  of  disturbance 
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correlations) ,  or  thst  the  mean  flow  with  sufficiently  snail 
disturbances  behaves  as  an  undisturbed  flow,  then  Q.  must  be 
independent  of  x.  (unless,  of  course,  l./L.  ~  0(1)).   With  this 
assumption  (diinensionally) 


Q.(l-,C)-  Q.  (xp  -  Q(X,Xj-,0 


(A.  6) 


and  the  mean  equations  become  (nondiraensionally) 


Continuity: 


^o-  La  c)  pU: 


Uo 


T. 


Ly  Po  '^^■  1    bpU 


'j  t" 


-    O 


(A. 7) 


Momentum: 


e 


^y^o^^x]  u.  ^.      _^i   (Jo:  Up,  ^x)   ^7,    '^^ 


U:  L: 


^X: 


"tr^/  ^  ^F, 


4 


^ 


2^ 


p  5,  L> 


^.' 


•X     Po 


L.  ^^^oV  ax, 


^ 


^ 


L^    ^ 


XL. 


L. 


-     /     .J 


/^^i^ij  ^y 


J 


(A. 8) 
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Energy ; 


ar 


'V 


c^.XAc:^^^ 


B  ^^)^v^^  a,  ^^ 


6/„  Z-,  J'  ^5c 


4- 


(^xY4^ 


.  /^ 


\o^) 


(A. 9) 


/:[M"l^^^^i^')^t^(^d'Si,2i^)]} 


At/:   ^ 


+    Jll  io   L^ 


1    3j        n  ^ 


vjhere  now 


A    -      ^^-^ 


^^    <5^. 


A    ^Q.l 


~J  u^    -I-  -J  U-,     .      -- 


a  ^  Lj    i  ax- 
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^ 


s 


^.-^Jax^       ^j 


4r   -^'  -i^   1 
J 


i, 


-for   -f  -  0(l) 
J 


From  these  equations,  several  interesting  observations  can 
be  Kiade. 

(1)  T\70  reference  mean  quantities  are  left  unspecified  by 
any  physical  or  geometric  consicierations:  v  ,  a  reference  normal 
velocity:  and  L  ,  a  characteristic  normal  length  scale.   The 
former  can  be  found  from  the  continuity  equation  by  requiring  that 
a  nontrivial  two-dimensional  solution  be  obtainable,  so 


The  latter  is  found  by  comparing  the  relative  magnitudes  of  terms 
in  either  the  momentum  or  energy  equations.   In  the  velocity 
boundary  layer,  viscous  and  inertial  terms  of  the  momentum  equations 
are  assumed  equal,  indicating  that 


L 


/ 


U.  L 


\Jn 


1     --      (?e 


7 


Alternatively,  equating  conduction  snd  convection  terms  of  the 
energy  equation  in  tlie  thermal  boundary  layer 


L. 


L 


/ 


Re^  Pr 
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Since  attention  herein  is  to  he  focuGcd  on  the  stability  of  a  flow 
field,  the  length,  scale  characterizing  the  velocity  boundary  layer 
is  chosen. 

(2)  For  the  assumption  that  the  mean  behavior  be  independent 
of  small  scale  variation,  Q  =  Q(X.),  to  be  consistent,  from  the 


moraentura  equations,  the  criterion 


tm'/}-' 


must  be  satisfied,  vjhere  the  amplitude  of  all  disturbances,  have 
been  assumed  equal.   If  1  represents  a  characteristic  wave 
length,  then 


Op 


<:<    0 


Re 


where  c<.   Re   >>  1.  Analytically,  the  amplitude  can  be  assumed 

y  '^y 

as  small  as  is  necessary  to  satisfy  this  criterion;  experimentally, 
control  of  the  amplitude  is  quite  a  different  story. 
Thus,  formally  assuming  that: 

(1)  i.ican  flow  is  steady; 

(2)  viscous  and  inertial  effects  are  equal  in  magnitude,  so 
that  attention  is  focused  on  the  behavior  of  a  velocity 
boundary  layer,  and 


U  .  Re    »  1 


Ly 
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(3)  all  disturbance  amplitudes  are  of  the  same  order  of 
magnitude  and  sufficiently  small  that  correlated 
disturbance  quantities  need  not  be  considered  in 
solvino  for  mean  variations,  or  specifically 


Qo 


Q. 


«    0 


(4)  flow  quantities  are  independent  of  the  z  coordinate 
or  the  body  is  "infinite"  in  the  z  direction,  so  that 


arj        B^  -  0 


(5)  (S°)-O0]and  (r^)  -  0(i)  but  WL,.   «  1,  the 
last  condition  implying  that  the  disturbance  wave  length 
is  small  compared  vjith  a  length  characterizing  changes 
in  mean  quantities; 

(6)  (1  /L„)  =0(1)  and  (1  /L  )"^  -  0(o(,); 

N    1  X   y  JjV 


y 

(7)  ^  =  0(1) 


(8)  dynamic  pressure  is  a  good  measure  of  the  system  pressure, 

so  p   ■--  n  U  ^ 
o    o  o 


the  m.ean  equations    to  0(Re      )    can  be  written   as 

Ly 

Continuity; 


-^ 


c)pV 


-  o 


(A. 10) 
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Momentum: 
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—  I  -     bu  -     A/y 


^'       sx     rSyy   ^y 


7  du 


(A. 11a) 


O    =      £>  P 

sy 


(A. lib) 


7^  I  n  o  w 


e  (^ 


4-    V    ^>'  1 


(A. lie) 


Energy; 


+  M 


•^^ 


^y 

z 

6yJ 


da  \  ^  Nl 


6y/ 


?T 


u 


<DX 


^  I  -^«= 


(A. 12) 
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Forraally,  the  disturbance  equations  are  obtained  as  the 
difference  between  the  instantaneous  and  mean  equations.   Per- 
forming such  an  operation  v/ith  equations  (2.1-2.3)  and  (A.3-A.5), 
the  complete  set  of  disturbance  equations  can  be  V7ritten  non- 
diraensionally  as 
Continuity: 

M  ^       ^^   ^^  ^»  ^  L,      ^.  u:   ^j  ^|.  °  J^       (A.  13) 


MoTTientiim: 


fi 


^  ^  -^^  r--    eo  ^ 


-  ^^^ ;  ^.J  ^ 


^„^/ 


?.  iy-  A  /'a  ^  ^  4-  ^  #  V  S' 


4-    2 


A  _:^    f^..    h.u    +    J.    a  ^H  ^   ^    hip   ~\+~]  lc¥ 


'^^cV,    i?, 


(A..U) 


i:)' 


7^  H  '\ 


J 


^i- 


''i'\ 
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Energy 


^^ 


rp 


~''  ^?e^"  ^>  ^^p. 


L^,    ^■]' 


£t         Vcp.  -^.J    /^     V   f^     ^  ^X 


^o.  T   i. 


^o^'  T.    J^-    ^  T      V    ^  J^. 


^oTYji^^/d/, 


^-'"^'^^S-^. 


ar 


-tx    y  ^' 


^_ — 


£i   ^4-5    - 


K'^  'o 


^i 


f/^o'  y  i;.  ix  y  i^ 


lo 


-,fz^;/;^^AA.V,r(S--A)) 


'X    L.>; 


-^f¥-^^). 


+A  Yii-^^ 


+    2/? 


(A. 15) 


/(^'j^'i',sr-/(^s^^'r>|)] 


X   ^z 


^^  +  i'  ^^  sfi  ,-,i  i;  [^  ^ 
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^■^'H'V   %      1j 


/       ^^-        CJ_I      I       f    '.1°    jiX       , 
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To  make  this  set  of  equations  tractable,  the  equations  are 
linearized  through  the  assumption  of  "infinitesimal"  disturbance 
amplitudes.   The  linearized  disturbance  equations  are. 
Continuity; 


^'4-  ^-  A  w  ^  +  ^^  ^'  il 


+ 


Uo    72  hi 


Ly    ^y 


(A. 16) 


e 


Momentum: 


-y 


f 


Bi. 


t/; 


dx 


uj; 


^■2-«^ 


4 


+ 


Ty  Uy     i;  hxJ     Lyisy     i.  5x/ jyJJ 


y 


(A, 17) 


k  1   ^^ 

Ll^.iJ     hi 


w„  i.  ^_y 


-   -^  ^"'   ■^"^^-  A  ^^  J 


-    1'  $ 


h  /oL^~y  ^  2>x       ^/ay       J,    sz^        L/'  ^^ 


^i;t)l^]^^l^f^^^^^ 
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^    aj?    ,    ^^  J   +   i/ 
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Wo 


c)  W      <5  - 


Ly    i^o     ^y 


7/    ^/ 


(A. 18) 
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_    r 


■^    \aJ  ^^    +  w,  4    ';  ^w 


dZ 


U,     L 


<^' 


7     -X- 


^  / 
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3/  V    ^y 


i'y   3; 


t 


-/ 


(A. 19) 


Energy: 

\  P  L  Wo)  ^t 
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hx   '    u^  A 
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Q>,To 


w.  i,   ,-^  £1  ^_   i. 


fT 
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^~> 
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Wo  ^.  w  ^ 
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■^   +-  A'   ^")  ^^ 


^y      Je  az/  ^4 
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Note  that  the  resulting  equation?  ^re  those  of  a  parallel  flov/. 

It  can  be  shown  bv  this  analysis  that  the  omitted  "non-parallel" 

-1/2 
terms  of  the  boundary  layer  equations  are  of  order  {<^      Re.  )    , 

Li  Lt 

y       y 

so   that   if  s^iall  values  of  c/       Re       were  expected,    it  would  be 

y       y 

necessary  Co  retain  these  teriiis,   Hov/ever,  since  it  is  expected 
that  heating  will  increase  stability  [1  -  7],  larger  values  of  this 
parameter  are  expected  than  for  the  unheated  case,  and  the  parallel 
flow  assu"ptions  will  be  even  ciore  accurate  for  this  boundai-}'  layer 
flow. 

From  this  analysis,  several  interesting  observations  may 
be  made: 

(1)  To  the  linear  order  considered,  there  is  no  provision 
to  account  for  thermal  fluctuations  of  specific  heat,  second 
viscosity  coefficient,  or  coefficient  of  thermal  expansion  (however 
slight  they  may  be) , 

(2)  Additional  length  scales  encountered  in  boundary  layer 

stability  can  be  easily  derived  from  dimensional  considerations. 

Consider  a  region  very  close  to  the  boundary  (viscous  sublayer)  with 

a  characteristic  velocity  u  ,  and  length  1  ,  such  that  u  ,  <<  u 

■^   ob        '-    yb  ob     o 

and  1  ,  <<  L   (all  other  conditions  remaining  as  previously 
yb     y 

specified) .   If  the  boundary  can  be  approached  so  closely  that 

%  -^X         27T 

>>      -  —— —  ^  then  equating  viscous  and  inertia  terms 


u  ,        L      c<, 
ob        y       ^v 
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of   the  raomentujTi  equatJ.on  requires   that 


7b         r^ 


I'lT 


Incidentally,  the  pressure  terras  is  also  of  the  saroe  orderj 
so  the  full  nomentum  equation  is  approximated  in  this  region  by 
Prandtl's  equation  [?-'/ ,p   6-9] 


2) u  ^  _  cV  +  ±   n^^ 
hi  dx       Re,   sy^ 


If  another  characteristic  velocity,  u.   .  and  length  1   are 

oc  yc 

"'         2tt 

chosen  so  that  1  << <<  — - —   and  still  1   <<  L  »  then  iu 

u         (X,  yc     y 

oc  Ly  jr  J 

this  region,  equating  the  inertia  and  viscous  terms  gives 


^]'  uj^ 


^L^PeLyl  iy/  Uo^J^ 


u  1 
o  y  c 

But  from  continuity  ', =  1  so  that 

■^   u   1 
oc  X 
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Vc      JTz^ 


Ly     V  c^uy^e^^, 
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This  region  corresponds  to  the  "critical  layer"  where  tiie  wave 

speed  and  mean  velocity  are  equal  (e.p,.,  u  =  c  for  a  two-dir,iensional 

K 


flow) 


For  example,  specifying  T  =  90"F,  T^  =  60°F  for  a  two- 


dimensional  f].at  plate  boundary  layer,  a  solution  to  equations 
(A.16-A.20)  with  homogeneous  boundary  conditions  is  found  to  be 
(see  Appendix  E) 


^^6 

=   7000 

-<R 

=    .122347 

^X 

=.-  0 

t^oD  =  3.68283  X  10~^ 

Cr   =  .210709 

C^   =0 

where  c„  =  u  at  y    =  .31852 
R        c 

Thus  for  the  length  and  velocity  scales  with  the  reference  length 
L  chosen  as  the  displacement  thickness,  6* 

u                             u 
°-»  51.355>>— '^^ =  4.7459  >>  1 


and 


u  ,  u 

ob  oc 


■/—     =   .085653 


1 

-2111-=   .19431 

6* 
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1  ^  "o 

In   the.  viscous   sublayer  0   <  y   <     — T-r~~  »    so   that  "^  > >. 

0"  u 

16.865  whereas  from  the  continuity  equation 


5 
U, 


2n  Re, 


,599.^5-7 


^      ^      -<Ls^. 


1  1 

V  C                 V  c 
In  the  critical  la}-*.--!-,  y  -  — ~-  <  y  <  y  -f  ~-  ■; ,  so 

u  26  26 

that  6.6978  >  >  3.7104  and  as  pointed  out  previously^  at  the 

u 

critical  point,   =  A. 7459. 

u 
c 

Thus  all  the  length  and  velocity  scales  derived  from 
dimensional  considerations  seem  to  be  in  accord  with  derived 
solutions.   Dote  that  these  two  length  scales,  representing  regions 
in  V7hich  viscosity  plays  a  dominant  role,  are  also  consistent  with 
those  derived  from  asymptotic  analyses  (see  Mack  [27,  pp.  4-4,4-5] 
for  example) . 


Appendix  E 

TherTP.odynaridc  and  Transport  Coefficients  and 
Tb.oir  Tenperaturo  Derivatives 

B . 1  Curve-Fitting  and  Derivative  Deternination 

Selection  of  the  curve-fitting  technique  used  in  this 
investigation  is  predicated  both  on  the  rcechod's  capability  of 
closely  approximatinp,  tabulated  data  and  for  predicting  high  quality 
derivatives,  with  particular  emphasis  on  the  latter.   In  the  ensuing 
discussion,  no  exliaustive  analysis  of  all  available  techniques  for 
achieving  this  objective  is  intended  or  attempted;  remarks  will  be 
limited  only  to  procedures  considered. 

Of  course,  a  given  set  of  data  can  easily  be  "fit"  as 
closely  as  desired  by  using  a  least-squares  polynomial  of  sufficient- 
ly high  degree.   However,  if  the  degree  is  too  large,  small  scale 
oscillations  about  the  primary  variation  'cill  be  evident  in  the 
resulting  curves  due  to  truncation  and  scatter  of  data  infomation. 
Utilization  of  these  curves  to  generate  derivatives,  then,  will 
naturally  indicate  even  more  conspicuous  but  completely  fictitious 
"wiggles"  in  the  computed  derivatives.   For  a  lower  degree  poly- 
nomial, not  only  may  data  be  poorly  approximated  but  variations  in 
higher  order  derivatives  night  not  even  be  seen.   For  both  the 
method  of  least  squares  and  other  such  empirical  formulations 
considered  (see  relationships  for  p  and  \.i   in  section  15.2,  for  example) 
it  is  significant  to  note  for  later  reference  that  property  varia- 
tion over  the  entire  temperature  range  is  described  by  a  single. 
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analytic,  continuous  function.   Thus  interpolation  is  not  betv.'een 
"exact"  tabulated  experimental  data  points,  but  between  approximate 
values  O.'hich  presumably  are  within  acceptable  experimental  error 
limits).   Derivatives  are  then  available  immediately  simply  by 
analytically  differentiating  the  indicated  equation. 

As  an  alternative  to  least-squares  polynomial  curve  fits,  the 
nethod  of  matherdatical  splines  can  partially  satisfy  the  previously 
established  selection  criteria  (assuming  availability  of  a  sufficient 
quantity  and  quality  of  aata  information).   In  general,  a  spline 
function  is  composed  of  piecewise  polynomial  arcs  of  degree  n  (where 
n  is  a  positive  integer,  n  >  3)  which  are  joined  at  prescribed  points 
such  that  derivatives  up  to  and  including  order  n-1  are  everywhei-e 
continuous.   Obviously  such  a  function  provides  continuity  of  the 
greatest  number  of  derivatives  consistent  with  the  iise  of  polynomials 
of  lov.er  degree  than  would  be  required  to  fit  all  data  points 
exactly  by  a  single  polyuonial.   Since  the  curvature  of  a  mathemati- 
cal spline  is  most  easily  controlled  v;hen  n  =  3,  or  a  cubic  spline, 
the  use  of  sinule  cubics  seems  to  be  an  attractive  v/ay  to  interpolate 
(in  the  strictest  sense)  experimental  data  points  representing  some 
physical  relationship.   As  applied  to  this  investigation,  an  explicit 
formula  for  this  "piecewise  cubic"  is  given  by 

J-l 
J-2 
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where  the  T.  and  the  J  (J  ^  3)  distinct  teroperatures  over  the 
normal  liquid  range  for  water 

0°C  =  T^<  T^  <  T3  <  .......<  Tj_^<  T^  -  100°C 

at  which  tabulated  property  data  (denoted  as  P.)  are  specified 
P   P  P 

J.      i-  »J 

and  where, 


0      -For  T^   Tj 
(T-T)'  ior    T>T, 


'J'    '  J 

The  J+2  coefficients  a,  b,  c,  A  ,  .  .  .,  A,.   are  then  determined  by 
requiring  s(T.)  "  P.,  j  -  1,.  .  .,  J  and  specifying  first  or  second 
derivatives  at  the  ends  of  the  tcnperature  interval,  T^  or  T.  [28], 
The  method  as  suggested  should  v/ork  perfectly  if  data  points  V7ere 
exact.   Not  only  vjould  this  data  be  exactly  satisfied  at  the  end- 
points  of  each  polynomial  segment,  but  the  methods  miniinization  of 
the  curvature  and  relaxation  of  the  overall  requirement  of  analyti- 
city  would  ensure  first  derivatives  of  high  quality  [29].   For 
reasons  intimated  above,  however,  it  is  unreasonable  to  expect  that 
data  lnrorrf:atii;n  cliould  be  satisfied  exactly.   Consequently,  the  pro- 
cedure outlined  above  is  modified  to  permit  departures  from  the 
prescribed  data  within  specified  tolerances  while  excluding  unwanted 
changes  in  curvature.   Numerically,  the  smoothing  is  accomplished  by 
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J  = 
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where  C  is  a  positive  proportionality  constant  and  r.  is  the  second 
difference  at  the  jth  point  defined  by 


J 


2 


T...  -^J- 


R    -P 


T 


+1 


P  -  P- 

J r 


Note  that  if  C  =  0,  then  only  the  first  sum  would  be  rainimized,  and 
the  procedure  then  corresponds  to  the  method  of  least  squares.   For 
large  C,  primarily  the  second  sum  would  be  minimized  in  an  attempt 
to  eliminate  undesired  fluctuations  [28],   Furtlier  note  that  in 
contrast  to  the  first  methods  discussed  in  V7hich  temperature  deriva- 
tives V7ere  obtained  after  smoothing  data,  this  one  computes  accurate 
first  derivatives  while  fitting  the  spline  function  to  the  specified 
data. 

For  the  present  application,  the  problem  is  thus  one  of 
mathematically  approximating  tabulated  fluid  property  data  suf- 
ficiently well  tliat  all  physically  realistic  variation  is  retained, 
but  spurious  fluctuations  are  smoothed  out.   Note  that  failure 


115 


to  acconplish  the  latter  should  certainly  be  revealed  by  evaluating 

higher  order  derivatives.   As  the  order  of  the  required  deriA'ative 

is  increased,  the  n\ore  difficult  this  problom  is  to  satisfy.   It  is 

just  this  difficulty  that  dcterrained  the  form  of  the  disturbance 

equation  adopted  for  this  analysis  (see  section  2.3.1). 

As  an  example  of  the  methods  discussed,  Figure  B.l  compares 

the  temperature  variation  of  the  first  derivative  of  c  obtained 

P 

from  various  sources  by  various  techniques.   Specifically,  included 
are  derivatives  obtained  both  analytically  from  least-square  poly- 
nomial fits  of  Touloukian,  et.al,  [30]  and  Kaups ,  et.al.  [7],  and 
numerically  using  a  cubic  spline  fit  with  smoothing  to  the  data  of 
Osborne,  et.al.  [31].   Consistent  v/ith  previous  discussions,  the 
latter  category  should  indicate  the  more  accurate  general  trend. 
However,  it  v/as  not  used  for  direct  input  to  the  numerical  program; 
the  following  discussion  explains  some  of  the  considerations  dictat- 
ing this  decision. 

(1)  Even  with  smoothing,  questionable  oscillations  in  deri- 
vatives are  indicated  -  specifically,  those  attributed  to  data  un- 
certainty.  In  short,  vrith  no  unusual  property-temperature  variation 
appearing  in  the  normal  liquid  range  of  water  (e.g.,  numerous  changes 
in  slope  and  curvature)  ,  data  can  be  v.'ell  approximated  by  a  least- 
square  polynomial  of  relatively  lo\7  degree. 

(2)  Since  the  aiaplitudc  of  a  fluctuating  fluid  property  at 
a  given  mean  temperature  is  assumed  directly  proportional  to  the 
first  temperature  derivative  there  (see  equ.  (2.35)),  it  is  essential 
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Fipure  B.l    First  tcnperatut p  dorivotivc  oi  c  o::   a  function  of 

P 
tciTiporaturc  as  obtained  U3ii;g  different  ptocedures. 
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to  know  v.'hether  and  where   erratic  behavior  in  this  derivative 
occurs.   Only  in  this  V7ay  can  logical  explanations  be  advanced  for 
aberrations  in  aigenf unctions ,  Reynolds  stresses,  dissipation,  etc. 

(3)  For  simplicity  of  evaluation  and  conservation  of  data 
storage,  it  is  nore  efficient  to  use  a  low  degree  least-square 
polynomial  or  empirical  carve  fit. 

Similar  analyses  v/ere  conducted  for  the  tliree  other  fluid 
properties,  p,  k,  and  \\.      The  results  suggested  the  property-- tem- 
perature relationships  for  use  in  this  stability  study  and  are 
discussed  in  section  B.2. 

B .  2  Property-Temperature  Variation  for  VJater 

Since  a  major  contention  of  this  study  is  that  a  fluid's 
property-temperature  variation  plays  an  important  part  in  establish- 
ing its  flow  stability  characteristics  in  a  heated  boundary  layer, 
accurate  values  for  thermodynamic  and  transport  properties  is  con- 
sidered essential.   The  literature  survey  prompted  by  a  search  for 
the  "best"  available  property  information  for  water  led  to  the 
following  results  and  conclusions. 

(1)  Rather  than  digressing  on  an  extensive  evaluation  of 
data  (v/hich  would  entail  an  analysis  of  experimental  techniques 
and  the  quality  of  information  available  from  each,  assessment  of 
systematic  error  limits  in  reported  data,  correlation  of  the  infor- 
mation -judj-ed  most  reliable,  etc.),  recourse  is  made  to  recent 
studies  of  such  a  nature.   Although,  as  evidenced  by  Figure  B.2,  by 
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far  the  largest  property  variation  with  temperature  is  by  the  vis- 
cosity, equal  care  is  exercii;ed  in  specifying  tb.e  behavior  of  each 
of  the  properties.  In  this  v.-ay,  discrepancies  betv^een  the  results 
obtained  in  this  numerical  study  and  those  occurring  in  nature  are 
less  likely  to  b2  attributed  to  inaccurate  property  information. 
It  should  be  noted  that 5  of  the  sources  encountered  in  the  litera- 
ture survey,  published  data  measurements  found  subsequent  to  t'nose 
used  herein  agree  to  vithin  experimental  error. 

(2)  Since  sons  of  the  "best"  data  available  are  generated  by 
polynomial  correlations  of  numerous  sets  of  data  (see  below  for 
c  and  k) ,  concern  was  experienced  regarding  the  capability  of  a  lov; 
degree  polynomial  to  be  able  to  predict  accurately  the  behavior  of 
higher  order  derivatives.   For  example,  v;hile  a  second  degree  poly- 
nomial equation  might  adequately  fit  a  set  of  data,  the  constant 
second  derivative  analytically  derived  from  it  may  not  be  realistic. 
Erratic  behavior  in  the  derivatives  can  also  be  expected  if  the  degree 
chosen  is  too  large.   Evidence  of  these  two  cases  has  already  been 
presented  in  Figure  B.l.   Using  the  smoothed  cubic  splii-^e  as  an 
approximate  standard  v^ith  the  data  set  of  Osborne,  ct.al,  [31]  (judged 

to  be  one  of  the  most  reliable  [30])  ,  and  knowing  that  only  the 

dc 
first  derivative,  —7=: ,  is  required  for  this  analysis,  one  can 

see  that  the  tliird  degree  polynomial  of  Touloukian,  et.al.[30] 

adequately  approximates  the  required  behavior  over  the  tcmpernture 

range  of  interest.   When  second  derivatives  are  required  as  v;ell 

(such  as  is  the  case  for  p,  p,  and  k) ,  the  above  procedure  regains 
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basically  the  same,  although  more  involved.   For  exaipple,  again 
using  as  an  indicator  the  sraoothed  cubic  spline  V7ith  data  recoramend- 
ed  by  Powell  [32]  (also  judged  to  be  among  the  most  reliable  avail- 
able for  k  [33]),  first  and  second  derivative  behavior  is  found  to 
be  best  ap-proximated  when  thermal  conductivity  is  defined  by  a  third 
degree  polynomial . 

In  short,  when  an  analysis  of  existing  data  uses  polynomial 
correlations  which  would  obviously  not  predict  behavior  of  higher 
derivatives  needed  herein,  information  characteristic  of  that  ob- 
tained from  the  polynomial  fit  is  used  instead. 

(3)  Regarding  the  required  property  variation  of  water,  the 
following  information  is  relevant; 

Specific  heat,  c  .   An  analysis  is  made  by  Touloukian,  et.al  [30] 
of  twenty  sets  of  data  for  the  isobaric  specific  heat  of  Jiquid 
water.   For  the  saturated  liquid,  recommended  values  are  computed 
from 


c  (c3l  g  ^R  ^)  =  2.13974  -  9.68137:aO  \  +   2.68536x10  \^ 
P 

-  2.42139xlO"^T-^       (T  in  °K)  (B.l) 


for  273°  <   T  <  410°K.   This  equation  is  found  to  fit  tlie  data 

analyzed  with  a  mean  deviation  of  .14%  and  a  maximum  deviation  of 

1.83%  (for  one  data  set  -  alJ  others  differ  by  no  m.ore  than  .5%). 

Of  particular  interest  for  this  study,  reference  [30]  indicates  that 

c  as  calculated  from  the  above  equation  "should  be  substantially 
P 

correct  to  within  one  percent"  below  400°K.   The  ability  of  this 


120 


cubic  to  predict:  first  derivative  behavior  has  already  been  dis- 
cussed. 

Density ,  p  .   Four  previous  works  are  evaluated  by  Gildseth,  et.al. 
[34]  and  coiaparcd  with  their  experiinental  results.   Recominended 
values  for  the  density  of  air  free  water  are  generated  for 
0°  <  T  <  80°C  by  a  modified  forra  of  the  Tilton  and  Taylor  equation 
(i.e.,  the  inclusion  of  the  exponential  term  which  is  negligible 
for  T  <40''C) 


r  /  1^   1    (T  -  3.9863)^(T  -}-  288.9414)     nii//n:    --  374.3  . 

(T  in  ^C)       (B.2) 


to  fit  experimental  data  with  a  mean  absolute  deviation  of 
.7x10   g/ml.   I.'hile  the  authors  do  specify  these  values  to  be  un- 
certain in  the  seventh  decimal  place  for  the  lower  temperatures, 
they  do  indicate  that  at  least  above  40"C  uncertainty  is  in  the 
sixth  place  and  increases  with  temperature.   Since  this  equation 
does  reproduce  the  densities  tabulated  in  the  International  Critical 
Tables  [35]  for  the  range  80°  <  T  <   100°C  (to  within  the  indicated 
uncertainty),  it  is  extrapolated  in  this  study  to  cover  the  full 
temperature  range,  0°  -  lOO^C. 

Dynamic  viscosity,  y.   The  experimental  measurements  of  ICorosi, 
et.  nl.  [36]  are  compared  with  ten  previous,  similar  sets  of  data. 
The  authors  propose  a  correlation 
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log  y  =  -  1.64779  +  ^  -'irTTs —  ^   ^'  ^"  ^"^   ^^^    ■*"  ^^''  "^^'^ 


to  fit  their  data  over  the  teraperatiirs  range  20°<  T  <150''C  with 
an  average  deviation  of  0.17%  and  a  maximum  deviation  of  O.A9%  at 
40°C.   However,  since  this  data  seeras  to  be  characteristic  of 
the  mean  of  the  other  coanared  data  for  a  specified  temperature,  an 
alternate  formulation  is  used.   Interpolated  data  is  generaged 
instead  by  using  (for  the  same  teiTtperature  range)  an  equation 
recommended  by  R.  E.  Manning  in  reviewing  the  paper; 

1    /  ^2Q°C  ,    1. 37023 (T  -  20)  +  8.36xlO~^(T  -  20)^ 

J-08  <.  ^    }  ^Q9  ^  ^  (.^.3) 

(T  in  "C  and  y  (20°C)  =  1.002  cp) 
This  equation  is  found  to  represent  Korosi's  experimental  data  to 
V7ithin  0.05%  average  deviation. 

Thermal  Conductivity,  k  .   Of  the  wore  than  sixty  sets  of  experi- 
mental data  available  on  the  thermal  conductivity  of  liquid  vater, 
Touloukian,  et.a]. .  [33]  evaluate  seven  as  being  the  most  relitible. 
They  then  correlate  these  and  generate  recommended  values  for 
273.16"  <  T  <  413.16°K  with  the  second  degree  polynonrLal 

lO^k  (ccsu)  =  -  1390.53  +  13.1937T  -  0.01903981^  (T  in  "K) 

This  equation  is  found  to  fit  the   data  considered  with  a  mean 
deviation  of  0.24%  and  a  maximum  deviation  of  0.82%.   As  indicatGd 
in  previous  discussion,  however,  a  second  degree  polynomial  appears 
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to  be  inadequate  to  describe  the  behavior  of  the  second  derivative 
in  the  temperature  range  considered.   Thus,  for  application  in  this 
stability  study,  a  third  degree,  least-square  polynomial  is  fit  to 
an  earlier  correlation  of  thermal  conductivity  data  by  Povrell  [32] 
with  the  following  result 

k  (pwatts  crf-'-K"-')  -  -  9.901090  +  .1001982T  -  1.373892xl0"S^ 

+  1.039570::10~''t-^      (T  in  "K)        (B.4) 

Values  generated  by  this  equation  deviate  from  the  preceding  one  by 
no  more  than  .25%. 

Due  to  the  nature  of  the  formulation  of  the  mean  and  dis- 
turbance equations,  It  v;as  found  convenient  to  calculate  the  pro- 
duct of  the  density  with  the  viscosity  and  thermal  condtu-:tivity , 
p  V   and  p  k,  rather  than  these  properties  individually  (see  equation 
2.20  for  example).   Thus  temperature  variations  of  c  ,  p,  p  p,  and 
p  k  as  well  as  their  temperature  derivatives  required  for  this 
analysis  are  shown  in  figures  B.l  through  B.5.   The  ratio  of  the 
second  to  the  first  viscosity  coefficients  is  computed  using 
equation  (3.3  of  reference  [37,  p.47]  and  the  data  of  Pink.crton  [14] 
and  is  included  for  reference  in  figure  B.6  (see  section  2.1). 

Since  the  stability  characteristics  are  found  to  be  sensitive 
to  assumed  property  variation  v.'ith  temperature  (see  Chapter  IV)  ,  a 
comparison  is  made  between  fluid  properties  calculated  by  equations 
(B.l  -  B.A)  (tjiose  used  for  this  jnvcstigation)  and  by  the  least- 
square  polynomial  fits  of  Kaups  and  Smith  [7]  (used  by  V.'azzan  et.al.) 
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Deviation  of  the  latter  from  the  former  (of  the  form     ih 

'-   "'^'  tcAUr'S^eJ.ol. 

(^)     X  /OO  »  for  example,  where  po  is  the  reference  quantity 
evaluated  at   T  =  BZ'F)  is  sho'ivm  in  Figure  B.7. 
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Figure  B.3   First  and  second  tenperaUurc  derivatives  of  p  as  a 
function  of  temperature. 
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Tigure  B.5  First  and  second  temperature  dcrivativp  of  pk  as  a 
function  oT  ter.pcraturc. 
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Figure  B.6  Second-to-first  viscosity  ratio  for  water  as  a  function 
of  temperature. 
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Figure  E.7  Deviaticn  between  the  property-toinpp.rature  variation  of 
Kaups ,  et.al  [7]  and  tliat  specified  in  Appendix  B.2. 


APPF.KDIX  C 
MEAN  FLOW  SOLUTION  VARIABLES 

For  reference  purposes,  the  mean  velocity,  temperature  and 
property  variation  through  the  boundary  layer,  as  well  as  their 
n  derivatives  conputed  bj'  the  program  are  plotted  in  this  appendix 
for  the  flows  considered  in  this  investigation. 

To  transform  this  information  from  the  n  coordinate  system 

to  a  more  physical  one,  Figure  C.l  has  been  included.   If,  for 

example,  the  reference  normal  scale  is  denoted  as  y   (i.e.  y  -  "-■  ), 

o 
then 


or 


In  this  formulation  the  advantages  and  disadvantages  of  using  n 
(rather  than  y)  as  the  independent  variable  may  clearly  be  seen. 

Advantages 

(].)   Using  the  mean  density  variation,  the  physical  step  size 
is  "automatically"  adjusted  within  regions  where  the  temperature  is 
significantly  different  from  its  free  stream  value.   This  means  that 
although  the  n  step  size  remains  fixed  throughout  the  entire  in- 
tegration range,  the  actual  distance  traversed  will  be  varied 
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constantly  within  the  thermal  boundary'  layer.   The  larger  the  wall- 
to-free  stream  temperature  difference  (and  so,  the  more  severe  the 
thcrriial  gradients  near  the  wall) ,  then  the  smaller  will  be  the 
.initial  steps  from  the  vail.   Admittedly,  liquid  v7atGr's  density 
variation  with  temperature  is  not  as  significant  as  miglit  be  desired 
for  such  a  formulation,  but  for  large  temperature  differences,  it 
is  still  sufficient  to  permit  selection  of  a  constant  step  size  out- 
side the  thermal  boundary  layer  while  adequately  describing  varia- 
tion within  it. 

(2)   If  information  obtained  from  solution  of  the  mean  flow 
equations  (2,20)  is  to  be  used  directly  in  defining  the  variable 
'  coefficients  for  the  disturbance  equations,  tlien  both  must  be 

expressed  in  terms  of  the  same  independent  variable.   If  the  latter 
is  to  be  solved  instead  with  y  as  the  independent  variable  (to 
compensate  for  the  first  disadvantage  discussed  below) ,  then  a 
constantly  varying  step 


1 


^'^itl  t 


must  be  used  within  the  thermal  boundary  layer.   It  v.'ould  seem  to  be 
inconsistent  (particularly  for  ]arge  temperature  differences)  to 
solve  the  mean  flow  equations  in  terms  of  the  llov;arth-Dorotnitr.yn 
variables  (eqn.  (2.14))  and  then  assume  the  m.ean  density  to  be 
constant  (p  B  1)  to  solve  the  disturbance  equations  (as  apparently 
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Wazzan  et.al.  [  3-5  ]  have).   Such  an   assumption  should  raise 
questions  concerning  the  accuracy  of  the  variable  coefficients  for 
the  latter,  especially  in  the  wall  region  where  density  changes  are 
■greatest. 

Disadvantages 

(1)  Outside  the  tliermal  boundary  layer,  the  n  step  size  is 

times  that  of  the  equivalent  "y"  one.   For  the  temperature 

range  considered  r,.,    is  alv/ays  greater  than  unity  (T  <  213^'F  with 

0*'  w  — 

T  =  60°r,  B  =  0  for  x]^.    <1  )  so  even  when  the  displacement  thickness 

CO  (^  JV 

is  chosen  as  the  reference  length  (y  =  6''0  >  a  larger  r\   step  would 
actually  be  required  than  for  a  comparable  y  integration,   llils 
means  to  obtain  the  same  accuracy,  more  steps  would  be  required  v;ith 
n  as  the  independent  variable. 

(2)  An  extra  calculation  is  required  to  transfoi-in  back 
from  the  q  to  y  coordinate  system. 

It  was  decided  that  the  advantages  outweighed  the  disadvan- 
tages and  so  the  problem  vjas  formulated  and  solved  using  q  as  the 
Independent  variable. 
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Figure  CI   Relationship  between  nand  "-^  for  various  wall  temperature! 
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Figure  C.2   Variation  of  n^^  with  changes  in  the  wnll  temperature. 
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Flp.ure  C.3     Second  velocity   derivative  evaltmted   at    tlio  vail   as 
a   function   of    the  wall    teir.pcrnture. 
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Fip.ure  C.^i   First  derivatives  of  velocil.y  ami   temperature 
evaluated  at  the  v;aJl  as  a  function  of  vail  temperature. 
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Figure  C.5   Specific  heat  boundary  loyer  profiles  for  various 
wall  temperatures. 
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Figure  C.6   Density  boundary  layer  profiles  for  variouf:;  wall 
temperatures . 
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Figure  C,7  Bounvi.nry  l.-'ver  profiles  of  the  first  derivative  of 
density  for  various  wail  temperature's. 
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Figure  C.8    Bonndnry  layer  profi'rs  of  the  Kecond  derivative  of 
density  for  various  wall  temperatures. 
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Fipurc   C.9      Boundary    layer   profilrs   of    tlie   donsity-viocosity   prnrluct 
for   varioiir?   v;.'ill    tenperatures . 
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Fl^,ure  C.IO   Boundary  layer  profiles  of  tho  first  clpriv.;il-.iv(?  of  the 
density-viscosity  product  for  various  wall  tcmpfjrature.s . 
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FijMjre   C.ll        Bountlat^/    iayor   profiles   of    thp    second   derivitive   of    the 
density-viscosity  product   for  v.irious  wall,    tcnpcratuics , 
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Figure  C.12   Eounrlary  layer  prol'llea  of  the  cien.si.t:y-thermal  conductivity 
product  for  various  v;all  temperatures , 
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Figure  C.i3        Boumlary    layer   profiles   oi    tho    first    (i(;rivatlve   of    the 

density-thermal    conductivity   product    for   various   v/alL    ten.pcratures , 
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Fif^ure   C.l-'t        liounclary    Jnycr   profiles   cf    Llie   fie^cond   dor ivauive   of 
the   dens J  t.y-thermal   conductivity   prodact    for  vnrious   w.iJl 
temperatures , 
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Figure  C.15   Boundary  layer  velocity  profiles  for  various  wol.l 
tenperatures. 
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Figure  C.16   Boundary  layer  profiles  of  the  first  derivative  of 
velocity  for  various  wall  teraneraturt^s. 
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Flfurp.  C,i7   Boundary  layfir  profiles  for  the  seccnd  dorivativc  of 
vcJocity  for  various  wail  temperatures. 
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Figure  C.18   Boundary  layer  temperature  profiles  for  various 
wall  temperatures. 
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Fir.ure  C.19   Boundary  layer  prcfilc;  of  the  first  derivative  of 
temperature  for  various  wall  tcnneratures . 
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Fip,uro   C.20        EoimcUiry   layer  profiles   of    tho  second   derivative   of 
tainperature    fo>:   various  v;nll    tonperatiiros . 


Appendix  D 
The  Computer  Prcp.ram 

D . 1  Prugram  Dss cription 

The  protrraiii  \-nritten  to  solve  the  stability  problem  forrrulated 
in  Chapter  II  consists  of  a  mapping  processor  (BLSMA.P)  ,  the  main 
program  (BLSTAE),  and  eleven  subroutines,  all  of  v.iiich  are  coded 
in  FOPJ'RAN  V.  A  brief  description  of  each,  along  with  its  role 
in  the  overall  conputational  scheme  and  references  to  numerical 
procedures  (where  applicable)  is  given  in  Table  D.l. 
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D .  2   I  np  j£t 

A  description  of  data  required  to  successfully  execute  the 
■program  described  previously  is  now  given  in  Table  D.2.   Further, 
the  required  forr.iat  for  this  data  is  illustrated  by  a  sarcple  data 
deck  in  Figure  D.l.   I.'ote  that  data  input  is  arranged  so  that  the 
first  card  lists  only  the  logical  control  variables  for  the 
program:  the  second,  all  rcean-flow  related  infornation:  and,  the 
renaining,  disturbance  £10;%'  and  stability  related  inforn-^ation. 
If  only  the  nean  flov7  solution  is  sought,  then  only  the  first  tv7o 
cards  are  required:  all  others  will  be  ignored  by  the  con^puter. 
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Program 

Car 
1 

d   Field 

Variable 

Description 

Propratn  Control  Parameters 

1-30 

KOTE(I) 

A  series  of  logical  variables  identi- 

fied in  BLSTAB  which  controls  the 

operation  and  output  of  the  program. 

Mean  Flow 

Parameters 

2 

1 

TWALL 

Wall  temperature,  T   (specified  in 

w 

degrees  Fahrenheit) . 

2 

2 

TINF 

Free  stream  tei-nperature,  T   (specified 
in  degrees  Fahrenheit) . 

du 

2 

3 

BETA 

Falkner-Skan  parameter,  M  - 

u  dx 
e 

2 

4 

EPSILN 

Specified  boundary  layer  edge  para- 
Bieter,  where  TT  =  1  -  EPSILN. 

2 

5 

H 

Step  size  for  mean  floX'7  equation 

calculations,  h  .   It  is  half  that  used 
m 

for  the  disturbance  equation  integra- 
tion, h,,  and  must  be  greater  than  or 

equal  to  n,.  /700  where  T]^    is  the 

0                  (J 

position  defined  as  the  boundary 
layer  edge. 

3 

Stability 

Pararp.sters 

Maxinuui  number  of  iterations  per- 

1 

ITIIHAX 

mitted  to  cornpute  an  eigenvalue. 

3 

2 

ANGLE 

An  orthonornal  basis  is  reformed  if 
the  angle  between  any  two  solution 
is  degraded  to  less  than  this 

quantity. 

TABLE  D.2   Description  of  Input  Data  for  Progrc-n  BLSTAB. 
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Ca-i 

-d   Fie]d 

Propram 
Variable 

Description 

3 

Stability 

Parameters 

(conto) 

3 

EPSLH 

Indicates  the  degree  of  refinement 
of  eigenvalue  estiinate.   Iteration 
is  terrdnated  when  F(ITER)  , 

P^(EV(ITER)_) 

RE(EV(ITLR+1))  '  ^"° 

AIMAG(EV(ITER)) 

AITlAG(EV(ITER-fl)) 

are  all  less  than  EPSLN. 

3 

A 

NEV 

Specifies  which  stability  parameters 
are  to  be  fixed  and  which  are  to  be 
varied  as  the  coiTiplex  eigenvalue, 
EV. 

^F       EV(ALFAR,  0>!EGA) 

2  EV(ALFAR,  ALFAI) 

3  EV(ALFAR,RE) 

4  EV (OMEGA,  RE) 

5  EV(OKiEGA,  ALFAI) 

6  EV(RE,  ALFAI) 

A 

1 

MM 

Kumber  of  eigenvalues  to  be  computed 
(maximum  of  10) 

5- 

■14   1 

REDLS(I) 

Reynolds  number,  Re.*. 

5- 

-14   2 

APJJLS(I) 

Real  part  of  the  wave  number,  a*6* 

K 

5- 

■14  3 

AIDES (I) 

Imaginary  part  of  the  wave  number, 
a*j.5*. 

5- 

-14  4 

FREO(I) 

Frequency,  w^  . 

All  of  the  four  preceedinf,  quantities 
are  the  best  estimate  of  the  1th 
eigenvalue  to  be  computed. 

TABLE  D.2      Continued 
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D . 3  Prog ram  Listing _ 

A  complete  listing  of  the  FORTRAN  V  source  deck  used  to 
generate  the  results  in  this  investigation  is  given  In  this  section, 
Note  that  whex"e  possible  j  the  routines  have  been  written  in  general 
so  that  application  to  systems  of  order  greater  than  six  is 
possible  inmiediately  v.'ith  perhaps  merely  a  ch.ange  in  the  DUiENSIOI! 
statements.   This  is  particularly  obvious  in  the  routines 
associated  with  boundary  layer  integration,  ADAI:IS  and  DIFFEQ,  and 
with  the  Gram-Schinidt  Orthonormalization  algorithn,  GSORTH, 
EGNFCN,  and  ORTHCR. 
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bLSTAB 

1  Di;.  C;,'j  I  Ci    C(1'J>701)  h'LP  i"<m  ,  7  O  i.  )  ,  Y(6.3'j!.  )  -LVCjO)  » 

^  lr<^L;L'-;(i.D)  .,-,t.-CLS{10)  .  AiuLSdO)  iFRfiaClO)  ?  ZQ  IbV  t  3  »  6  >  .V- U  r  !M  ^5  >  :5  »  351 )  f 

3  ^;f-"i.LiCTG(3.o5n  »(-U;jCC0»351) 

'»  '           OuUBLL    PaCCISICIJ    HMIA 

5  COMPLEX    r,LVfFLUCTG'"OIST»P 

6  U^..GICaL    iJuTECfHP.OVAK 

7  EvJU.lVALE'iC£    ( P  .  FLljCTO )  f  {  ZOIST  f  PIJ; 'C  ) 

a  CO.Vf.OiJ    /:^^Kl/ALPt^Ari,ETA^p■|<IlJr■Vt;LI<2/R^r,.'5I^JF('l•)  »C??/DLK3/G(tt)  » 

9  iGF  ['■•')>  yjpp  (  4  )  /uLK^i/LLPAR  »  C/f.;LK5/K'  r  f-,i  KTOT  .  !<MAX  >  KMAXPVBLk?/ ALF  AR  r 

1 0  2Ai.F  A  1 1  Kl  >  Gi>;GA/ L■■U^o/^:on:  (  30  ) 

11  HcAD     lb,  200)     (i,OT£(I)  I  I  =  lf30) 

12  u.-ne   i6ri.i6)    ti^ofEC  I) .  r-i»30) 

13  116    fjh'.-i\'\     (IX»3UL1) 

I'i  2U'J    FoKMaT    ioul.U 

lb  C  ir-    lOsICAL    VAKlAf-Lt;?,    -.OTEd)    -    TRUE,     ThE    FOLLOwIMG    ACTIOil    IS    TAKEM5 

16  C  (iOTEC..)          TEHi-ilMATE    PKOGIvAM    EXECUTION    AFTEl    CALCULATING    AigO    STOiaNG 

1/  C  li-i    UMU    a    'LAM    l-LOv;    PAf'A.VclERS    AMO    STAtjiLlTY    E&N    COEFFS. 

i.S  C  l,0TEl2)          CAuCULA'lE    OijLY    EIGEliVALUE,     If    FALSE,     EIGENFUi-XT  lOi.S    ALSO 

19  C  CO'IPUTED. 

20  C  |;0TE(3)          i'lU  TE    ETA » VELOCITY,  TEMPERATURE,     AND    THEIR    EIA    DERIVATIVE?. 

21  C  liCTECr)          V.RITE    KTOT    Ar:0    ETAuEL 

22  C  i.OlU'))          l-^vITt    C(I,K)'S    -    STABILITY    EOUATlOl-j    COEFFIC  lEK'TS 
2.i  C  i,OTE(t)          ■.■,,'aTE    ALC-LEf  ITR!.;AX,EPSLrJ 

ii^  C  I.OTECn          i.'i'^ITE    ETADLS    AiJO    RATIO    OF    H.L.    TO    DISPL.    THICKr^ESS     - 

2d  C  ,\0-|£(!j)          .■;R1U:    EIGEr;V/LUES 

20  C  ;^01E(V)          -UTE    "LAN    FLO.-;    ,<OUi'JOAKY    LAYER    STEP    5I£E  ,  TjALL  r  Tli  ;F  ,  ETC  . 

27  C  IvJILiiO)       ..kITc    GUESSES    OF    I-'EAfJ    TEi'P    AIIO    vELOC    Gr;AD    rtT    V.ALL    AMD 

2.1  C  T.v.Ii'    I'.CRE^Et;TS;     ,\RITE    A'-Vr:-PT0TIC    ERROR    I|.!    SOLUTION,    E. 

29  C  .,OTE(U)        flEAO    .'•■EAtJ    FuOv.'    A,-!D    STABILITY    EGli    COEFFS     IliFORr- ATTOrj    FROM 

30  C  STut<AGt    U.-ilT    (J.     JF    FALSE,    P.v.OGRAM    KECOhlPUTES    ThIS    INFO 

31  C  pruCR    TO    PFvOCEELINo    TO    EIOEiiVALUE    Pri03LEM. 

3.?  C  ,jOTE(.12)       C0,.:'>'JTE    A^.'0    V.RITE    RELATIONSHIP    L-ET;;EEN    ETA    AMD 

a  c  y/(uIsplacen;::it  THici';i\'tss) 

3't  C  ,.'0lEll3)        :.><X^'-:    SPECJFIC    LT  ,  uENS  IT  Y  ,  I  GENS  I  TY  cV  ISCOSI  TY  )  » 

3b  C  t..'E:jSUY*ThEI  :''AL    CO.'jOUCT  I V ITY )     A':D    IliETK    Ef  i\    DERIVAiIVfS. 

3l>  C  ,,01EU4)        .-..urE    ITEHfEVdTER)  »GA:-;-.A,S['-.!:a,Z(J.- ItKPTOI  )  'S 

37  C  NuTEU;^;       CO    PUfE    LENSITY, /ISCOSlTY,A-iL'    TIlLRiiAL    COl.'D    FLUt TUATXOr.S 

oO  C  ;.ijTE()6)       C.;,.rUTE    KEY..'OLi)S    STkESS 

3'-J  C  i.OIcCiV)       .■.•JlTE    UEK;-    -    -.1.'.    OF    CMMC.-.OR::aL  T /"?  A  T  i  0,/S    RE'JUIRLD. 

4i;  c  i,o;e(Io)     ;.-<irE  x(j)'S  =  Cw:-'„inI:  c  cor.FFiciENfs 

m  C  .lOfEli'))       .'IITE    1(l,i)'S    -    CIGE:  rUr;CTlo;  S    VALUE    AT    WALL 

4<i  L  i;DlE(2':')       .■.'.•irE    F(lT-r)    Ai.i;    CAIIS  (F  (  I  f;- !0  ) 

'43  C  ,.OTE(i:l)   '  :,iCATc.  T.,/.T  IT  ;<''.AX  E,VCEEU.-:n  AND  PEST  EST  OF  EIGri, VALUE 

'+'♦  c  !.'jte(.?.2)     r-   mr'.Ax  kxce'/(..':i),  write,  all  ev(Ite:<)  and  f(ite;-;) 

••t)  C  NUTE(^.3)        -.'UfL    ^;.R^:;^^)     -    KOSITiO.iS    OF    O'iTno.  JONImALIZAT  IONS 

^(i  L  i,0Tt.(2'H        US^    P.'.iJPERlY    V ',(■  t  ATION    G(I)-(,(U)     OF:     TOULOUK  I  Aii ,  ET  .  AL  ,  J 

M7  C  GIL^CcTHtET.AL.  Jr.OnOSI  ,E  r.  Al,  i     A'  D    PO.VLLL,     IF    fAL'.E,    USE 

nii  c  I, 'AT     or   K.urs  c;  S;;ITm    (fi.-..:;   IF   i,orEC5(i)    F,..LSr) 

4'J  C  ,,uTE(._!.0       .■.•■a.E    CAt-,'.(r  riT  .■)  )     =    VALi.'L    AN'.j    POblTU.N    CF    ELST    SCLUIICN 

!JU  C  ,.0-|Ll,:..)        I,;<,LU::E    I.E.   S  11  C    EL'  CTUA  i  IC.  j'',    ill    SOLUTiOIJ 

bl  C  N0Ii-(27)        J;n-!C,.IE    UZ-i.    uE    rULLElVS    IT^|;^TI^N    SCHEME 

S2  C  I,u1E(2o)        !   !C_L:E    TE--i'i:.''«M","E    FLUC  TUAT  )  lir:  ■<    IM    SOLUTI'jr;     (ASSl.l'.'ING 

bi  (.  NJ'lt.(3u)  =  Tr,'J:. )  .     I?    FALl^E,    1:0    TL'".P    FLUC    ARE    A'.'.SU  •,E0 . 
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5=+  C      i,0T£{2'J)       HKlur    AlsbLr:    fiETV.CEii    SOLfJ    VlCTORS    on    first    iTERAriCi) 

bb  C      .iOV£.(iU)       iJb,:    LO'„lLL«S    l";L'aTI0;JS.     if    FALSEp    use    f.OUATIGfiS    OF    'JAZZAMf 

bo            L  0<,,.-,U;vAt     A.'.D    SMI  fH, 

b/  Mr    (J 

bd  N    =    3 

69  K..-,AX    =    Vul 

bO  R£AD     (b.J'jU     TV.ALL/TI.ir.BLlAfCPSlLUfH 

6i  ijl    r.Jr..''iAr     '. -^r  \  ij  ,2  il.l(, ,  d>  0  !.n  .  c> 

■od  C    K-  A  A    =    y.AX    Mo.    OF    STtPS    ALl.O./ED 

b3  C    11    >    C,.;    -    t:T,.JLL/(^^■.AX--l)  ■  Oi<    Olr'lEMSICIJ    STtiTS    fiUST    ijE    CliAMGEO    FOR    Ky.AX 

o'f  T.iALL    -    i...>lL    +    'if-?. 0.9 

ub  Tii.F    ::    T  L  if    +    'tLyy.b9 

6b  AlPHA    -     (l,£TA    +    l.Li)/2,0 

67  R    =    TliK/l'Jl  ,69 

6U  S    -     (T. .,-,'.. U"ri:jF)/US).l. 69 

6y  C^    =    (T„ALL-TIf,F).fb./9. 

yu  Call  tli-v;ai<(o.OiO) 

71  Ir     (^.0TL(9))     ..'RlTE    (&'1(;9)     H»  LETA  >  rv'ALL-'t59.59f  TluF-(|b9.  60  f  PRlNFf 

72  11,-EFS1L.-J 

73  119    FoiOiaT    ('    STEP    S12£    -'(iUO.a/'    FALHsKiER-SKAM    PARArViETER    ='fF10.';/ 
7'4  1<    v,ALL    T;;i;rERATURE    -'»FjO,4/<     FREE    STREAM    TEI-'.FERATURE    -',FiO,H/ 

7b  2«    FREE    SlREAVi    PRA;.ijTL    IJO.    -'»Fin„fa/»    GOUNDAt<Y    LAYER    EDGE    DETERMTIJE 

76  3D    ijY    U    ='  i.r  11.7/) 

77  I,-:    (hOTE(li)     .AnO,     .not.    NOTE(I))     GO    TO    1 
76  u    C„LL    Lr-EO.UII»EPblLr)»TvvALLrTIIJFfETADLS) 

79  ETACEL    -    iM.PAKdrKTOT) 

riU  Ir     (L'OrF(..;>)     vRITE     (f..-10.1)      CKf  (GLI'^R(I  fK)  »  1-1:7)  fK  =  lfKTOT) 

61  icl    FuRHAT    (/b  IX'  '  »  '  f  15X.-  '  f  »  '  »37Xf  •  '  '  '  15X;  •  r  f  '/3Xj  'K'  f  1^X«  'ET/i'  ,i6Xf 

6Z  ;.',•.'  f  17a  (•.,'  ,  1'jXi  'vv'  r  17X>  '  ThETA'  <  iax»  'TH^TA'  »liiX.  'THETA'/ 

b'^  IF   (;;oT;_{.f))    aritf    {Gf1o2)   ktot-etaci:l 

6b  l.j2    F.-KIiAl     (/•     i-TOT    ='.Ii+/'     c'lAOEL    -'»E13,3/) 

6o  jtf-     (i.Gr::.»7))     .■.;<iTr.     (6, loo)     CfADLS 

67  l-.o    F;;r.'-;M     (/'    tVA    SU;;.:;CRIfT    uELTA    STAR    -''Elo.S/) 

b.5  .■RiTE     '-c)    rRTi.;r  fETAfl.SfET,u'JEL»KT0T;'-;0TE{2i*)  f  I !.  T.'ALL  •  T  IIJF  .  KETA  » 

o9  iF:,-SlLI-'  (  <Ni  •■'A.;(I  /i;)  .  f-l  ,7)  .  1', -1 .  KTwT  )  ,  (  (  C  ( 1  -  K  )  f  I-l  f  I'O  »K-^lfKTOT) 

"ju  IK    (:,'0T':(i))    STOP 

91  GO    TO    2 

92  i    R^.sU     (i-i'":-!-)-.^.)     FRIIiF»r.TAOLS;ETA0i:L7KT01  iH20VAf-;>HA,TWALLA>TIi:FAfnET 

93  lA;.JtF5A.  I  ;  ;LPAf;(  I  .  r, )  »  i  -  i  r  7  )  >  K-1  i , ,  I  r.  r  )  r  (  (  C  (  I  ;  K. )  f  i -1  ?  I'U  »  K-1  r  K  fO  T  ) 
9^1  IF    (iior.i.Jt)     .Ai.b.    .■!2JV:,n)     GO    TO    J.  0 

yb  i-    {•::!F(2;)    .op.  h2o^/^i'.)   go  to  ii 

96  -.0    Ir     (iiiih;  (  C^'--m)  f)  .005)     .£0.    0     .-'..•-■j,     T'vlVt  C  (TV.ALLA-TWALDh  CTIi/.-A" 

97  lTi,,F)  t  {..r  iA  .-DLrA) +  (r::'S.,-cPSU.I.')  )  .  ).  ,E6)     .EO.     l.i )     GtJ    TO    9 
96  _!  I    ,;  ;IT::    U.;  117) 

99  117  r.,,,.,AT  ;/!:<.'., o'.  •+•' J/bX.  '  Tm;-:  REGUE.-.rru  I)r0R.-.;ATU);j  IS  t;CT  THE  SA^'E 

100  lA:,  Tii.j  SI. .;■■_;;  I'^i  ui;n  <  . '/ix^''.i.i(  >  -  •  )/> 

Itil  STi;f^ 

102  'J    i;.-  (;:;j;r.(b,)  IRJT,-.  (nrlji)  (  K  f  ( -"-LI  ,.:•■  (  I  »  K  >  »  1  =  I  »  7  )  r  K-1  f  K  TO  1 ) 

103  IP  (.juTKl'.)  )  /.RIT'.  (i,>;v,2)  KTOT  » V.  1 /i':EL 
1U+  IF  (•.■jL:(/))  ,.!alr.  u-»;.h:>)  ETAOLS 

I'Jb  2  l;-;.'>!,i  r^'l-J:-.)     X1R:'>AAf  A'^'.l. ''.»LPSLI,'ir.JV 

loo  2.ri  fJ^tr\^     (';  UM-'IO. -t  >i.] 'J  .•:  /  UG) 

107  IF  (.juCtl?)  )  y,!<H(.  <^flf'3)  i  K  .  IC  (  )  .  :■ )  .  1 -1 1  I'O  .  K.;:  1 1  KTOT  r  b) 
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108  loo   F-'OlsMAT    ( (  r+f  i'*t;o.3) ) 

109  I.--    (,;^rf(oLj)  )    MilTE    (6»ion 

i.io  I:--   (,;,0T.   rjGTL;(oO))   \,i'-;rrE(6(ioa) 

ill  ].;••     (..N'JT.     f,0TQ:(<i3n     .ViU  TE  (  6  r  1 10  ) 

Xl^  lu7    P0I-n:aT     (/<     t.+  t    KcSULTS    FO.'i    rrUATIOr.'S    of    LO.VcILL    +*t') 

11^  iOii    FuKMAT    (/'     **^■    f.ESULTS    FOR    EaUATI0;j5    Or    V.'AZZAiM'    Ol<AN;UKAr     M\0    S--:ITh 

H'^  1    jf*^'/) 

lib  liO    Fv.Ki-iAT     (>     ^i.^    V.ITHjUT    1  E.M.-^l"ATU:7C    FI-UCTUAT  ICjIIS    *(■*'/) 

lib  IF     (;J0TE(o))     ivKITF    (6il0'l)     AfiCLE »  I  TRMAX  .  EPSLN 

llV  .iu4    Foiu-iAT     (/•     I'.IH    AliGLE    AL.LO„EO    nETV.'LSri    S.OLM    VECTORS    :::'.•  FIO  ,  4/ 

lly  i'    NAX    i-.'O,     iTEKATXCijS    PlK'N'ITTFG    ='rl3/ 

iiy  2'   co:jv'epol;jce  ciUTLraon  For?  eigeuvalues  ='»eio.4) 

120  M;,lcP  =  -i;.0*ii 

1^1  K:-!AXa  =  (K,-;AX/2)+l 

li^''     C  EbI"I:-'.,-,TE  i!E»0^;iiGAf  A!-iO  ALFA  To  USE  IN  E.V.  CALCULATlO;i 

laJ)  REAO  ('o,2Gj)  :-'i;-»  (HECl-Sd)  ,ARDLS(I  )  .AIDLSd)  »FKE0(1)  fl^l.MM) 

12t  2ji    FoKMAT     (I-t/(';Ei6.y)  ) 

12b  Do    3    l  =  lM-i,M 

126  i\'jTE(2l)    =    .Tl^UE. 

127  RE    =    REDLSm/ETAriLS 

128  AuFAI    =    AIuLS(I)/ETADLS 

129  0..-,EGA    =    FREO(I)*RE 

130  ALFAK  =  AKOLSCD/ETAOLS 

131  IF-  (MEV  ,E&.  1)  EV(l)  =  C:lPLX(ALFARf OMEGA > 

132  IF  (iJiiV  .EO.  2)  EV(l)  =  CMPLX  (  hLFAR  » ALFAl ) 

133  IF  <;.iEV  .EO.  3)  EV  1 1 )  =  CMPLX  ( ALFARf  RE) 

134  IF    l\iE'J    .£0.    'D    EVll)    =    Cr-iPLX  ( Or-.EGA  f  RE ) 

135  IF  INEV  .EQ.  5)  EV(l)  --.  Ci-'.PLX  ( OMEGAr  ALFAl ) 
13b  IF  (,'<EV  ,E0.  6)  EVd)  -  C.-;r-L.X(l?E»ALFAI) 
137  Ev(2)  =  (1.01.U.G)*EV(1) 

130  E7(3)  =  CXPl.Xd  .0(;Li*REAL(EV(.n  )  »1.02*AMAG(EV(1)  )  ) 

139  Call    ulS  FES  Vdy^i  TP,.,AX  » l  iSTEP  f  Ai  'GLE  f  Y  >  EPSl.i; .  PR  INF  .  ME  V  ?  53  »  P  r  ZOIST  ) 

140  FrEO(I)     -    O.v.EGA/RE 

141  AruLS(I)    -    ALFAK*i.TAOLS 
1^+2  H^OL^il)     ~    is^E*t.TA.-.l.S 

143  AIOL'o(I)     ::    AEF  AI  i-r  I  ADLS 

144  C;  -    Oi-EGAvALr,  K/(A!.FAf.-r.*<:-^ALrAr^*2) 
14^  C;  -  -G'EGA-'  "LFA'/(  A>.-FA '^ffZ- ALFAI*^-2) 

146  IF  (,;OTt(.,))  ..l^;V-  !od03)  AR0L5  ( I )  »  A  IDE?  t  U  '  kEOLC  {  t  )  f  n;EO  ( I ) 

147  liCRfCl 

14«  lutj    F>.,R!w",T     (/•     ALFA.-.    -'.Elu.e;/'     ALFAI    -•»tl.;,.3/'     RE    ::'»E10.8/'     FReO       = 

14y  I'rclu.o/'     CR    ^'fEl^-.'/'     C.i    -•■El«.o) 

lio  IF   (  .jiE(.d)    .OR.    .;.0T.   ■;0Tr.(^l))   co  to  3 

Ibl  Ki'TOT    '-"    .'.lOT/^  +  l 

1^2  Ii-     (,;,ur.    ,.,OTE(.lS)     .C'^..     .NOT.    ;,OTl.(30).)     00    TO    5 

153  kV.xITt  (Or  1.  j.c::) 

Ib'^  lj,2    F■.^l^^^Ar     (//I       KP'  .'-X.  '[lE'.Sn  Y    hU'.'CTUATI  0., '/ 7:<  r  '  VISCOSITY    FLUCTUATJo 

iS'J  l,i      TpK.RI-'AL    C0i;U(iCTr:lTY    FlLCTOATIO;;      ArjS'CE'.FLU)          Ai;.o(VISFl.  U)          A 

lii5  2BS(TCr-LU)  '  ) 

l'^7  Ov-,    4    .N^I.KPIOT 

Ibd  Ti;.'-;P    -    i•L.^■u^(5>.:V:.-'l) 

Ib'J  C„lL     rE';VA,'{TEJ  P»1  ) 

IbO  I|.     (., Ill', .(,_'■,,)  )     FL'JC  fO  (.•!/.■-)     "    O(E')  tGI'(2) -Y  (3.K) 

161  Ir     l.il^'T,    ;:OTE(?.A;;  )     FLOC10(2fK)     -     (  0  .  d  r  0  .  C  ) 
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.lo2  Fi..UCTG(3f;;)    -    G  t  3  )  kGP  (  3  )  ♦  Y  (  3  » K  ) /G  (  ?. ) 

163  Fi_U0Ti3i'W,N)     ■•:    G  (  4  )  *bp  ( 't  )  *Y  (  5  fK  ) /i",  (  ?.) 

io'i  li-     (iuT  (T.,.\LL.-TIi\f)     .r.^L.    0)    V.'iilTE    (6*111)    K  .  (  FLUCTG  (  J  »  K  )  ,  J-?.;  4  )  f 

165  1  (^Ar>S(FLUcr.O(  J/a)  )  .  J~2)H) 

loD  ill    F0iM':A7     (  m.Mi'l'i.O/aXfL'.E.l'K.'J) 

l'j7  H    COinTIIjI-.'e: 

loa  b    If     (,iiOT.    ,-iOTF.  (16))     GO    TO    3 

io9  C       CO;  i'.JU:    T£^.;S    OF    £.,,;:;:,C-Y    G>>LAmCE    EGUATION    FOR    TiEUlKAL    DlSTURc  ANCt.S . 

170  C         Ft',\\:ii.fK>     ~    .\f-.Y,;OLOS    STRtSS    PRODUCTION 

171  C         f-L..i^(<':fK)     "    FJ'.OLUCTlor.    OVZ    10    VISCOSITY    FlUC  TiJAT  I0(  .'S 

17c;    c   Func(3.K)  =  i'.-:t:sr.uRE  emcuoy  productioh 

173     C    Fu.-iCCtjK)  =  jIS5IPATI0;t 

ly'i     C    Fu.l;:  (5,-K)  :;  jISSlHATIO;i  DiJE  TO  t,or:COriSTA:!T  ;:EAN  At.O  FLUCTUATIOfl  OE.":, 

17b     0    Fu;i^(6rK)  =  .'I^OOUCTION  OUL  TO  ^'EA^i  VISCOSITY  GRADIEmT 

176  Iv.alE    'w-.n^O 

177  ii4  F^Ki'AT  (/'  E'-JEhGY  P;;ODUCTIOfJ  ANO  Dlf^'.ilPAT  ION  TF.PyS  FOk  ticUTRAl.  STa 
I7u  IbiLITY  (AlFAI-0,0)  •/'  KP  RE  STftESS  P:^00  VI5C  FLUCT  PROD  PRES 
179  2:S    PROD' .DAT 'OISSIPATIOf;' »oVp 'DEM    DiSSIP            VISC    GKAD    PROO'/) 

l;iO  Do    6   K,=lrKPTOT 

Ibl  L    -    2*K-1 

Ic'^  FUKC(1»K)     =    ~ULF"'AR(3.L)  !■  ALFARt-RF-i-rsi-;  AL  (  Y  (  2  ,  K  )  "i-COf  i  JG  (  Y  (  1  (  K  )  )  )/C(l'L) 

lc3  IF     (.NOT.    ■.)Tt(50))    FUI.C  ( 1  .•  K  )    =    FU.lC  ( 1 .  :<  )  *C  (1 .  L  ) 

lb4  FU:nC(2m<)    -    ~L:LPAi:(3fL).J  ( i?EAL  (F1..UCTG  (3  »  K  )  ^CONJG  (  Y  (  5  »  K)  )  ) /C  ( 1 ;  L  '  *t^a 

liib  i+AEFAi<**^  +  AVr-AG(ColvJ6(Y(l.K)  )  *FLUCTG  ( 3  M<  )  )  ) 

IGo  IF    (..\0T,    ,i0fc(3D))    FUi;C(a»K)    =0.0 

1«7  Fuf;C(3'K)    -    -ALFAI-*RE*FEAL(COIUG(Y('VrK)  )  =i' (  (  0  ,  .•  1  .  )  :;<FLUCTG  (  2»  K)  ^ 

y bfi  i  (  kLP Aii  5  2  » I. )  ~C[: )  ^  C  (  1 .  L )  ■<  C  <  1 ';  t  L )  /2  .  * Y  (1  ( .K )  )  ) 

l.j'J  IF     C.tiOT.    i-;OrE(30))    y-UiX'.O'l'.)    =    0.0 

100  Fu;iC(4>K)    -    -C  (2f  L)f(  AI.FAi-;»'f'>  +  Y(lf|<)-!CO'lJG(Y{i»K)  ) '' C  (  1  r  l.)  ^i-e^- 

191  ICOI.JO  (  Y  ( 5 ,  ,0  )  '^  i  (  ^j  f  lO  /C  (1 '  L  )  ^^R-^2  .  ■i=/,l.FAR**2*  A  IkAG  (  COI!  JG  (  Y  (  fj  f  K )  )  * 

VJ^  2Y(i»K))) 

1V3  Fu;iC(fjfK)    =    J'..*C(2(L)*C(7»L)*ELPAR(GrL)/C(UlrL)*ALFAR**2*(  (  ( 

Vj  ^  lli^PA;;  ( ?-.  •  U )  -CF, )  -i  C  (  i  »  L  )  )  •■  -i;  "FLUCTG  ( <'. »  K  )  •■  COf '  JG  (  FLUCTG  (  2  r  ,\  )  )  l-  (  C  (  1  'i  »  L  )  / 

19'j  i:::  .  )  *  <  2*  Y  (  1  r  K  )  *C0;JJG  (  Y  (  l  ,:•.))  +  (  BLPAr^  (2(E)  -CH  )  >:--C  ( 1 :,  »  L)  -i  C  (  I  »  L  )  + 

1-Jo  .1 A  I  f-AG  ( (.0M..:6  ( FLUCTO  ( 2  .  K  )  )  -'  r  ( 1 »  K  )  )  ) 

19/  IF     (.FOf.     .,0T^(33))     FU:c(twi;)     =    U.O 

I':  3  ForC  (cwl',  1    -    -■i.-»ALF.Ai-:  **--'  tC  (  2  r  I. ) /C  ( 1 »  L)  •+  (  C  (8  »  L  ) +C  (  i't  f  L  ) /2  .  )*Ali:AG( 

V)  <  iCC);<>JG(Y(l,lJ  )  vY  (2,K)  ) 

2UU  IP     (CF;     ,Gi,    tEPA'v(2»L)  )     00    TO    7 

2i!l  If     (C.%'Gr  .,4-PA.'(2,L-2)  )     ;  ;U  TE  (  6  . 1  lb )    K-1  i  K  r  BLf^•^R  (1 '  l.-H  )  (oLPAfi  (  1  ^  L  ) 

202  U;>    Foir'AI"     (•     i>;E    CRITICAL    POIiiT    IS    PE^E    V-HEKc  '  <  I^  f  •  <K<  '  >  1 3 1  •     OR     '» 

203  IF  (.2,  •  <:eT.,<'  'I  IF.2) 

204  f    .,  iHE     UwU^;     ,-'.)  {!  lj;,C(  J/K)  .  J-lf6) 
2'j;.i  ij3    FvivFAT     (  I'Jf  i^^Uij.'W 

207  ::■    Co.,Tll,OE 

•7.i.''J  ■      El, J 


169 


MFEQN 

1  SUbROUIINE  MFEQN(H»E.P5ILN.TWALLrTINF»ETADLS) 

2  D.i.MEiJSlOfJ    y  (15)  »OYUb)  fETMEND('+)  »TEST(t)  »Bl7f701)rC(m,701) 

3  1»GN(4)  »GI>iN('+) 

f  DoUbLE  PKECISIOU  Y f DY » H r ETA . X r Z » UELX » DELZ » SUMODD» SUMEVN 

5  LOGICAL  NOTE 

.6  CONiN'.ON  /riLKl/ALPHA»l3ETArPRIMF/BLK3/G('t)  rGPCl)  »  GPP  ( '+ ) /BL  K'+/B»  C/ 

7  lBLKb/i'J»K»KTOT»KMAX»KMAX2/DLK8/NOTE(30) 

8  EXTERNAL  OIFF 

9  DATA  I4U/2/,  (ETAENOd)  « 1  =  1 »  2  ) /2  .  » 10  .  /  .  (  TEST  ( I)  »  1  =  1 »  2  ) /I .  r  1  .£"8/ 

10  1»ETEST/1.0E-12/ 

11  X  =  l.uuu 

12  Z    -    -l.ODO 
li>  1  =  1 

l"*  MARK  =  1 

lb  K  =  1 

Ifa  i  ETA  =  O.ODO 

17  Yll)  =  O.ODO 

18  Y(2)  =  O.ODO 

19  Y(3)  =  X 

20  Y(i*)  =  l.ODO 

21  YCb)  =  2 

22  YC6)  =  O.ODO 

23  Y17)  =  O.ODO 
21  Y(d)  =  l.ODO 
2b  r(9)  =  O.ODO 
2b  Y(10)  =  0.000 
27  Y(ll)  =  O.ODO 
2a  Y(12)=  O.ODO 

29  Y(13)  =  O.ODO 

30  Y{m)  =  O.ODO 

31  Y(lb)  =  l.ODO 

32  Call  AUAMS(15»H»ETA.0.YfDr.l»DIFF) 
ii)  GO  TO  (i+f  30)  f  MARK 

3t  30  B(1»K)  =  ETA 

3b  B(2»K)  =  Y(2) 

^ii  E(3.K)=Y(3) 

37  BC+rK)  =  DY(3) 

38  B(5»K)  =  YCI) 

39  ■   B(6fK)  =  Y(5) 
"+0  B(7»K)  =  DY(5) 

HI     C  COMPUTE  MEAN  PROPERTY  VARIATION   THROUGH  THE  BOUtJDARY  LAYER 

42  Do  1  1=2.4 

43  G.jd)  =  GP(I)*B(6»K) 

44  1  GiJN(I)  =  GPP(I)*IB16.K)**2)  +  GP(I)*B(7rK) 

4b  IF  (ri0lE(13))  WRITE  <6f.l01)  K  r  G  ( 1 )  »  (G  (  I )  »  ON  (  I )  *6(I )  »  GK,M  ;  I )  *G(  I )  r  1  = 

46  12r4) 

47  101  Fuf<MAT  (I4.10E12.7) 

40  C  COMFUTL  MEAN  FlO'J  COEFFICIENTS  FOR  DISTURBANCE  EQUATIONS 

4y     C  FOK  MlCH  PKESbUHE»StT  \/iS2=~2/3.  FOR  THERMO  PRESSURE.  SET  VIS2  =2.1 

bU  V152  =  -2./3. 

bl  DVIS2  =  0.0 

b2  C(l.K)  =  1.0/G(2) 

b3  C(2»K)  =  &(3)/0(2) 
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bl  Ct3rK)  =  {2.+VIS2)*(2.*GN(2)**2-GfMN(2)-GN(2)+GN(3) )-GN(2)*DVIS2 

5b  CC+rK)  =  GN(2)*(2.+VIS2)-eN(3) 

bb  C(b»K)  =  -GP(3) 

57  CCbrK)  =  12.+VI52)*(GP(2)*(3.*GN(2)-GN(3) )-GPp{2)*B(6rK) )-GP(2)* 

bd  1L)^/IS2 

by  C(7rK)  =  -(2.+VIS2)*GP(2) 

6U  C(8»K)  =  -GN(3) 

6i  C(9»K)  =  -GPP(3)*B(6fK)+B(3.K)~GP(3)*B(iWK) 

62  C(lOfK)  =  -6PC2) *(1.+VIS2) 

63  C(ll.K)  =  PKItvlF*GCl)/G(t) 

64  C(12»K)  ;:  -GNIU'+) 

6b  C(13pK)  =  -2,U*6N('+) 

6b  0(14, K)  =  2.»GtJ(2) 

67  IF  (ll.-V(2)) .LE.EPSILN  .AWD.  Y  C+J .LE.EPSILN)  GO  TO  13 

6a  6  K  =  K+1 

69  "+  CALL  AL)AMS(lb»H»ETA»l»YfDY»l»DIFf-') 

7U  5  GO  TO  t6f30) >    MARK 

7i  6  IFCETA.LT.ETAENDCI) )  GO  TO  t 

72  7  Bll  =  r(7)*+2  +  Y(9)+*2  +  Y(8)+*2  +  YC10)**2 

73  Bi2  =  Y{7)*Y(12)  +  Y(9)*Y(1'+)  +  Y(8)*Y(13)  +  Y(10)*Y(15) 
.7"+  B2i  -  bl2 

7b  622  =  r(12)**2  +  Y(lt)**2  +  Y(13)**2  +  Y(15)**2 

76  Cl  =  ~l (Y{2)-1.UD0)+Y(7)  +  Y('+)*Y19)  +  Y(3)*Yi8)  +  Y(b!*Y(10)) 

77  C2  =  ••'(  IY(2)-1.0D0)*Y(12)  +  Y(45*Y(li+)  +  Y(3)*Y(13)  +  Y(5)+Y(15>) 
76  UEM  =  Bll*fci22  -  021*812 

79  DELX  -    (C1*B22  -  C2+B12)/DEM 

80  DELZ  =  (C2*B11  -  C1*B21)/DEM 

81  IF  (|40TE(10))  WRITe(6.10'+)  ETAEND  ( I )  r  XOELX»  ZtDELZ 

82  lot  FORMAT  (FIO .or tD20 .8) 

83  E  =  (Y(2)-1.0aO)**2  +  Y('+)**2  +  Y(3)**2  +  Y(5)**2 
8'*  X  =  X  +  DELX 

8b  Z  =  Z  +  DELZ 

86  9  IF(ABS(OELX/X) .GT.TEST(I) .OR.ABS(DELZ/Z) .GT.TESTd) )  GO  TO  o 

87  10  IF(E.Lr.ETEST)  GO  TO  12 
66  11  IF  (I  .EQ.  NO)  STOP 

69  1=1+1 

90  Go  TO  3 

91  12  MARK  =  MARK  +  1 

92  TEND  =  Yd) 

93  ETASTP  =  ETA 

91+  IF  ((JOIEUO))  V.'RITE  (6»106)  ETEST»ErX.Z 

9b  106  FORMAT  (/•  MAXIMUM  ALLOWABLE  ERROR  IN  ASYMPTOTIC  SOLUTION  ='»El3.6 

96  l/«  ERiiOR  IN  FINAL  ASYMPTOTIC  SOLUTION  ='»ElJ.6/«  FliJAL  ESTIMATE  OF 

97  2:    VJALL  VELOCITY  GRADIENT  = • » Dlb.8/22X f • WALL  TEMPERATURE  GRADIENT 
96  3  -' fDlb.B/) 

99  IF  (N01E(13))  ftRITE  (6rl07) 

lOU  lu7  F0Rf-iAT(/37X'  '»'fliX.'..'  »23X»  '  r  '  r  1 1  X  >'.»'»  22X  r  '»';  lOX  ,',»' /2X  >  'K  '  r 

101  13Xf  'Sf'LC.Hr.  '  »5X.  'DENSITY'  ,tiY.,  'DENSITY'  ,5X»  'DENSITY'  rbXr  '  (DEN*VIS) 

102  2'  »2Xf  '  (GEN*VIS)  '  »3Xr  '  (L;Er>i*VIS)  »  .bXr  '  (DEN*  TO  '  fOX.  '  (DEN*TC)  *  »3Xf 
1U3  3'(bEN*IC)') 

104  GO  TO  3 

lOb     C  ENSURE  KTOl  IS  ODD  FOR  STABILITY  EON  IMTEGRATION  AfJO  ETADLS  SOLUTION 

10b  13  IF  <.K/^*^    .EQ.  K)  GO  TO  8 

107  KTOT  =  K 
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iUd  C    C,.,LCU^ATl    ET,i    SUliSCKli'T    DlLTA    STAr!    =    CTADLG    USI;;G    CO.'-'PO'ol  Tt£    S.TMPSONS 

109  c  Kv,ur:  fok:;lL.', 

110  If     (  i;il  (T..,AUL-"riMF)     .Mi-..     0)     GO    TO    I't 
Hi  DcJiiNi     -    !'.1A 

11^  ETAQLS   =    CVASrf-    -   YEr;o 

iij  Gj   ro  i? 

11^  14    5j-',0:)U    -    .1.0 

115  Sji.t'ViJ    =    0.0 

116  I.si:    r    KTOT-? 

117  I.il    =    Kfyr-l 
116  D.v    1":.    K-i,  i:^?.,2 

119  lb    Sj,-.Oi>l    -    S'J.-.OUui    +    Cdf'O 

120  Dw    16    Kz2r  [;-il»£ 

121  lb    So'r.El'.'M    =    Su:-;i£Vt-;    +    Cd.K) 

l^;^  Ut-liI^JT    =    iCdfi)  fCCl<i'.TCT)i-2.:.^SUyOL;D>'+.-*SU(.;EVii)+ETA/3./lMl 

125  eiACLS    -■    ijc;f;IHT~Yl-.l:D  +  £Tfi,SlT~ETA 

ISi*  1?    Ir'    (fiOTEC/))     '..'iilTE     (6rlO?.)     1 . -EPS  .'.Lr;  j  DENTi  IT/ET  AoLS 

12b  lo2    FuKI'iAT       (/'     CLOU'JCAKY    LAYcK    TtilCKiJESS    CJ-"- '  «  F6 , 7  .  '  ) /DISPLACEHEl:T    Tm 

12b  llClM'JESS)    -'rEJa.O/) 

127  IF  (.,'*0T.  i-jOTLHi;'))  RETURN 

12U  KrurN'.2  ::  KIOT  ~  2 

129  Yi,o;;o:i  =  0.0 

iJO  W,aT£     (6rlU3)     YhOi  DM  »  B  ( 1  >  1  ) 

1^1  luJ    FoKr:ar     ('     Y/dTISPLACEMtNT    THICKTIFSG)  •  f  6X  '  '  ETA  ' /2E20  .  8) 

1^2  DO  la  K=ifi<T0T;;2'2 

li3  Y.<0..0-'i    -    rriOMDi ,+  ( ij.  (1  ,  K-i  2)  -b  (1  f  K  )  )  *  C  C  ( 1  r  K  )  +'i .  'KC  ( 1 »  K  +  l )  'rC  ( 1 .  K-i-2 )  )  /6 . 

1^'f  l/ElAULS 

155  16    WKITE    (GfiO'S)     Y,\0;  :DM.  B  C 1  »K+2) 

IJfi  1C5    FuK?;.\T    (2E20,0) 

137  RLUJaii 

15U  Eii'O 
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Uli-F 

3.            C  Sj,UOuTi!;£    TO    CaLCUI.ATc    TriL    UIFFtfx-l'.'iTI AL    EOUATIOilS 

i:  SUi.KCjn.ic    ClFF(i;TA»r 'DY) 

3  D^-USLf:    PiuiCISIolj    ElAiY.OY 

"*  Dir.Lf-.biO:^;    Y  (lb)  fOYd'^)  rGXm)  >  02  < '• )  '.-fMS)  »6iJX(2)  ?Gri2(2) 

.5  CoK:f-iO;<    /.'■i.M/ALPMA-fjrTAfPKINF/tiLI'-J/GCl  )  ,Gpl').)  rC-PPCl) 

6  Tc.i'.P    -    Y(0 

7  C.^LL    TL- ■>,';,;;(  rENpil) 

8  0.)    1    I--l,4 

9  G,-.  (I)     -    bT'd)  *Y  (9)  ^ 

10  ■      1    G/.(I  )    -    ,;p(  I)  ,:Y(1<)) 

11  Do    2    I-J?(i 

12  G/(I)     =    G:^;T)-Tir(5) 

13  G,.XCI)     :-    i;-r-MI  )i-Y(b)  -YC'.')  tGPO+YdO) 
IH  2    Gi;:(:)     ::    .rPF(I)*Y(L))  .'YC).ii)+GP(I  )  :  Y(15) 
16  D((i)     =    Yli^) 

16  Df(2)     -    Y(o) 

17  OUi)     '-■    ~•o.i(^)*Y(5)-;l>^TA*(1.0/6(2)-Y{2)*^2)^■/,LPHA  +  Y(l).+  YC3)  )/G(3) 

18  Dr('4)    ^    Y(5) 

19  Or(L))     =    -0,-.(i^>)*Ylb)"P;ai;Ft=G(l)+ALPHA*Y(l)'!-YC5)/G('>) 

20  0( (6)     =    Y(7) 

21  0Y(7)    -    Y())) 

22  OHQ)    =    -■..X(.5)*lyY(5)~GliX(3)*Y(3!-C.i;(3)*Y(6)  +  (bl::TA*(GX(2)/6(2) 

23  iH-^.0  +  Yt2)-li  (7))-ALPHA+(Y(6)tY(3)  'rY  U  )  *Y  (8)  )  ) /G  (  3) 
ZH  DY(9)     -    Y(J.O) 

25  0  .'  do  )     -    -g:-:  !  4  )  *DY  ( ■"• )  -6:iX  ( 4  )  *Y  (  5 )  -on  ( t|  )  >Y  ( ]  (l )  --PRltjF*ALPHA*G  (1 )  * 

26  1  CGX  d  )  *  r  ( J. )  "  Y  (  5 )  -i-Y  ( b  )  *Y  (  5  )  ^-Y  d  )  ^  Y  d  0  )  )  /G  « '^ ) 

27  Didl)    -    i'd2) 
20  DY(12)    =    Yd3) 

29  0  f  ( 13  )    =    -G/!  {  3 )  -HOY  C  3  )  -Cf  iZ  (  3  )  *Y  ( 3  )  -or,'  (  3  )  'Y  (13 )  !-  ( [ih"TA+  (  GZ  (  2 )  /G  ( 2  ) 

30  l  +  ^i.O*Y(;:)  r-;  (U/)  )-ALPhAv{YCl],)*Y(3)  tY(l).iYd3)  )  )/G(3) 

31  DY(1'^)    :-■    i(i.r,) 

o^  DY  d'-> )  =  -u-i  ( '^ )  +DY  (  V  )  -G:;2  (  4  )  +  Y  ( 3 )  -C!  !  ( 't  )  tY  d5  )  -PR  I  ,'iF*AlPI  !A*G  ( 1 )  * 

33  l(G2d)  fYd)^Y(5)+Y(ll)*l'Cb)+Yd)*Yd5))/G(4) 

3'+  Rcl  URN 

i'j  End 
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TEMVAK 
1 
•  2 

4 

b 

6 

7 

6 

9 

10 

11 

12 

1-i 

14 

15 

16 

17 

1« 

19 

2U 

21 

22 

23 

24 

25 

26 

27 

2a 

29 
3U 
31 
32 
33 
34 
3b 
36 
37 

3a 

39 
40 
41 
42 
4  3 
44 
45 
4  b 
4  7 
46 
4  9 
50 
51 
52 
5J 


SaoIvOUT  Xi;,-.    te:;vaR  ( "I  €IV  ,  VAOf\) 

L,,G1CAI-    I'iOTf-: 

Cy:..Mo.i    /■iL^2/KfS'Si::F(4)  » C2/BLi;3/G  C4 )  #GP(4)  i  CPP  (4) /OLKa/fiOTL"  (  aO  ) 
1  /:,;_l'.  1  /.'■LP.  1,1 «  Lf:  f  A  r  P;.,  Ii\F 

Qii'.'-;/:,IO  ;    m(21)  ..3(10) 

D„1  A  (A(l.  )  ,L-J.  ><il)/l  ,4C33Da9,0.."0  7^:501  fO.;:5?,a9602»0.00392,'W0.4(  159j 
11(0.  2oo']y  /;; ,  0  .  0.53 'i  o;)9  ,  45  .  2  jgavoM  f  137  .  19-'53uo3  f  14?.  .  7970  /'4  -i  i  b9  .  819762 
23,9,90-toi,Vif2.5a20i^'6i6.70562dlf  3,  7:SU4471?  0.61634 17  »  7  J,  376906.200. 
j7  t74'.s:'ic  ,  l'J7  .7604o/'b<6b  .c"6261i>bi  7,477945ri/ 


D...T/ 


(L  (L)  .L-1  .  19)/::.  lZy7'fr--9,6f;137E-o.2.i.)Ol53ut-.-5,-2.421,39[-: 


13  ,  9bo3  '  :^-.-o  ,  9414  »  5oo  ;'29  .  i?  f  6H  .  120&3  r  Ci .  0)  1445  ,  J7'l  ,  3  '  1  .  00,i , 
.<a.o6;_-H,  10'.'.  '  ,4342y44af-9,901G9.  .  1  fjC  193A »  "1 . '573fi92!:-4  »  1 
IF    1N„LX    =    0    C0-iPUTl5    r'.Lr.    STHEAM    VALUES.    GiIJF    (KAUPS    fi    Sh\l 
CU;-.P,,TF-S    OI.;F/G''.eF    Willi    P-iJF    .^UAi  IT  IT  li-l^    EVALUATED    AT    T    -    32 
IP    Ii'L-cX    >    (.1    CO-  PUT,-:    Piior-EPTY    VARIATION.    G(l).    AT    Sl-'ECIFIEO 
G(J)     -    Sf'ECJFIC    PEAT    AT    COr  STa;:T    P[<ESSUI^E 
GC2)    u    PZ,;SirY 
G(3)     =    i:.[:,.bi  VVtVISCOSITY 
G(4;     -    DPl^l5I  rr»(Tri£;-v;-;AL    COnCUC  riVITY) 
ALL    Tr.ESE    OUA,.TITI£S    APE    fJCr.'OIM.    W.R.T.    REFEREl-jCE    VALUES    AT 

T    =    S*TMi-.p+R 

Tc    -    (  r-l)=r-273.16111 

TPv    -■    TC    +    273. 1& 

If-     (TO    .Li',    l.U)    TEXP    -    0.0 

IF  (TC  ,GE.  1,0)  TEaP  =  EXP { -p ( 10 ) /TC ) /TC**4 

IF     (  .Pf-'T,     NOTE  (24)  )     GO    TO    1 

G(l/    ^    r- (1)+B(2)  tTK-;-i.i(3)  *Tr<  +  *2t-B('l  )*TK**3 

G(^)     :;    J  .."CrC-lJ(5)  )*  +  2+(TC4B(f.)  ) /li  ( 7 )  /  (  TC  +  H  ( 'i )  )+B(9)-:'TF. 

G(3)    -    0(i)  tj(li)/10.  »*(  (d(12)  •<(TC-20.)-tB(13)=.'(TC-20,)^^ 
IBCl'i  ))  ) 

G(4)     -    0(2)  +  (B(16)-Hi(17)  ^TK-^B(lB)'^'TK.**2-t-B(19)*TK**3) 

G.^    To    3 

1  G(l)     :;    A(  U -A(2)*T-HA(3)*T**2 

Gi.;)    -    A(',)-!-A(5)*T-A  (b)  *T-f»2  !-A(7)  i-T  +  *3 

Gi^)     -    1  .,j/(.\(c.)-,- (P)  ■  T  +  A(10)*T**2~A(ll)*T-*-  +  3t-A(12)*T*  + 
G(4)     ~    ",\(i3)->-A(14)  ^-T"A(  15)  f•T  +  ♦2^  A(16)*T**3 
3    II      (  U.JE<    .GT.     0)    GO    TO    4 
D,j    2    1  =  1.4 

2  &lr.P(  i'    -    S(  I) 
Ir     (. •iC'I  ::(..:.•)  )     P,U:i-    -"    M  .  B4  *G  (  1 )  ^  G  (  3  ) /G  (  4  ) 


1  .  "-57023. 
.03957E-7/ 
TH  FORM 

F) 

TEMP 


[IIJF 


XP>TC**4 
;<2)/(TC  + 


4) 


IF     (  ,i.'-'f  .    ,.0Tl(24)  )     PPI 
lA(c:'))'>'t^^'  3  +  .,(21  )  )'T*"4) 
RETUP.r, 
C .-.  ;PU ,  E    ( CG  (  i )  /fj. ) )  /G  ( I )    Apr; 
(  !...iL..-'l  :;:■-  ) 

4    ir    (,;;CT.    :,;irL  (;■»)) 

G,-  (  i  )     -     (..,(-<)  -I  ,;r  .-.    I 

C,J(2)    -    (-(TC-i- (L-)  ; 

1  rC'v,;  >  (TC  i'  2^  ■•.  (lO;  .. 

t'vi-  {  C)      ■■      ;   (   (-U  .  V'  I.  H 

io;i.) )  •  r^4,:(  ;i.ij^-(-  ) 
GP!3)   -■  vi'i;  )-(.(;. 

1C^/(:J(1'..)   i    (>     (14)   HC)   ) 


13. GG/  (  A  ( 17  )  "A  ( IC  )  *  r-i  A  ( 19  )  *T  +  *2" 


(U2GC1)/I:M2)/G(I)     ivPEI-.E    W    -    (T-TIf.'P)/ 


GO    TO 

(.:.    re 


B(';)  iTh.  i*;M+C2/g(1  ) 

rciP.  ;p(G)-•u^)  )/P(V)  n. 

(  vf  !  ;  )  i  C2'  ■■■(2)/(TCHUr-)  ) 
,  :-.  )  •■;■'  .  <  ■;:  ( '..•  ))/•=(  7  )  -i  !j  (  '! )  *  I-  (  1  0  )  -f-  (  ( :' 
■')/?,(::)  +CP'-*2-2.  •  G!-ia)  :^.r')/(TL-t 

c'do)  '(VL--?o. ) 


G(2)i!.U'J)  *TFXP* 
(10)  •■?..* 

i(.')) ) 

A.IO('-Mli)  iG(a)/u(3)  )  )'r 
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136 
57 

bci 

by 
60 
6i 

ta 

65 

6-t 
6L> 
66 
o7 
bO 
69 
70 


6i^p(5)     •-    uPi- CSJ+GPCM  +i'£~CP(^)>5>i2-2.  +  (0(lJ)  ^C2*  +  cV0  t  15  )  •)  (GP(3)- 

Gp  en    =    OP  (  ?.  )  +  (  J  (  17  )  ^2  >  rn  ( .U) )  x  TK  1-^^ .  ■"'  3  (  10  )  *"l  i<*-t  2  )  ^  C2  +  G  (  ?. )  /G  (  H  ) 

GPr'C'O     =    :,PP{i!)  K'3,*GP{a)  *  (CP(-t)-GP(2)  )•!(?..  ^\]{XS)  +6  , +0  (  19  ) +TK  !  .|-C2 
l*Ti+-G(a)/G('>) 

GO    To    6 
)    G,'(l)    =    (-rt(2)+A{:>) '.?..='T)-5/G(!. ) 

GpC?)    -■     (  .N(b)-2.  t.^(6)  -  1 +3  .  :■  A  (  7 )  -t-T +=i'2 ) -SS/G  (2 ) 

Gr-P{2)    -     (-2.  fA(6)+6<  ^A(V)  *T)*5  +  +  2/0(2) 

o,-^(3)    -    -G(::.)'i'(-A(9)+2  =  ^(^A(iO)  *-T-3.  +  A  C 1 1 ) +T**2+'+ .  *  A  (12  )  +T  +  *3)+S 

Gi'P  (  3)  =    2  .  *GP  (o )  +*2-G  ( 5  )  *  (  2  .  *A  (1 0  )  -6  .  *A  (  U  )  *  r+12 .  *  A  ( 12 )  i-T**2 )  >S-i^  +  a 

GPCJ)     -     (A(14)-2.*A(  lb)*T!-3.*Atl6)  *Tt.>2)*S/G(U) 

Gi'Ptil)    -     l-2.*A(15)+6.*A(i6)*T)*St-  +  2/G(a) 
)    D'.j    1     I-".  ;'4 
'    G(i)     --    G(  n/GIiil-  (I) 

RcTURil 

END 
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ADA,.,S 


1  SUuROUTiNc    ADaV-.S(!;»HjX,  I'JCTr  Y-OY.  Ii;pf:X.F) 

■  ^  D.;l.t.LL    P;'.£CISjUM    P  »  Y  t  W,  »  OYlI.  r  uYLLL  ,  I'T  L  .  OY  ;  DYR  ,  C;.' .  C-^  »  CH  ,  X  i  H 

<i  Dli-.rlK'blO')    -(  (lb)  .DV(.Ui)  fl' (I'j)  .  YR(i;.)  >D^LLL(l'j)  -triLl.  (15)  »DYL(15)  'OYfi; 

t  KlL.)  .C^(i;;)  <c:5(15)  fC^KlD) 

5  1    Ih ( JScT.Gf  ,0)       GO    TO    6 

■6  ^    I;- (liiL.ex.Lt.O)       GO    TO    ^ 

f  o    i\    "    (^ 

8  GO    To    ^ 

9  't    K    r.    1 

10  ;■)  c.LL  r-(x,Y»nY) 

il  Kt.'lurUj 

1 ^  b    Go    Y 0     ( 1 0  »  7  »  C . 9  J  11  )  .  K 

1J>  7    K    =    3 

i4  Gj    To    10 

ib  6    K    =    'i 

16  60    TO    10 

17  y    K    =    i3 

10  iU    Z-j    lOOl     I"1»N 

19  lOul    P(j)     -    Yd)     >     (H/2.0IjO)*OY(I) 
Zo  CALL    F  (X-iri/2.0(jO»PfC2) 

21  Oj    1002    I-lfti 

^l  100-2    Hd)     -    Y(n     +     (tl/2.0O0)=?C2(I) 

<^3  Ca.L    F  (A+;i/2,01.'0rP»C3) 

2it  0^.    lOU-i    TirlfN 

20  iOu.-i    P(  i  )     -    Y(  !)     +    H4C3(1) 

26  Cvi.i.  r- ( x+t;»r'f c't) 

27  Do    lOu'4    I-J  .il 

2a  lOo'l    Y.v(l)     •■•    Yd)     '.-     (M/o.O,Xi)<  (uYd)+2.0D0':C2d)+2,0D0.-iC.m)+C'KI)  ) 

29  C,>LL    F(X  +  ,lfYP,  rbYR) 
iO  Go    10    12 

31  11    0.)    llOl    l-\i\\ 

iZ  1101    lid     -    Yd)     +     (n/24  ,DDO)^  (!;)5.0D0*DY(  I)-rj9.0D0trjYL(I)-!-37,0O0i-OYLL(  .t 

6i  l)-9.-;;ji'vDYLLL  (  d  ) 

j'+  C;\LL    F  (X  i-ll(i-'>OYK) 

y^  Oj    .1  )02     l-l  ,M 

30  11(.,2    Y:<Cd     ■-    Y(J)     +     (ll/2'+.0L'O)*('?.0u0>OYR(  1)  fi9.000*nY(I)"b.ODO  +  DYLd) 
jl  d    CYLI.d)) 

3a  C„Ll.    r(.\  (-liiYli'OYR) 

:v->'  1^    X    r.    /l-M 

't.j  Do  Koi    i-i»;j 

'H  Y(l)     ::    Y.Ud 

M2  Dn-LL(d     -    CYLL(I) 

't3  D  1 1.1.  (I)     =    UYUd) 

4'*  ufl.(  d     =    OYd) 

^+5  12  jl    ,;i().  )     -     y/Kd) 

h6  H^.TU;■'l^ 

'+7  E.iO 
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Ul'..TtiO 

1 

•  2 
3 

b 

6 

7 
8 
9 
10 
1.1 
1?. 
IJ> 
)^ 
15 
16 
17 
16 
19 
2U 
21 

25 
24 
2b 
26 
27 
26 
29 

:ii 

v52 
i6 
34 
35 
36 
37 
36 
39 
'10 
'H 

42 
•i'o 

4u 
'1 7 
4y 
49 
bO 
bl 
b^ 
b3 


K.'-'Ax » H ,  ANGLE  f  Y . K^SLt; » PK Tr  F ( t ;e:v ,  i 
; VI A I  s  1 6^;  ^  ?  n  nt  >  o  y  i ;  it  i  z  v  f 


1^\0 
21 


SiJbP.O'ul  I  Uc.    0 1  ST  t  (3  ( t;v ,  I T 

Co^  ;PLt.X    Z ,  At;- ;\  » i> .  Y  r  X ,  G,'i 
if^E!-0;-.iC.lroc.»H2(D0 

L;..GICAU    il>jTt-:».v.Ui.LLK 

CO,.,iVOf-.    /tn.:\4/oLPAK  r  C/r!LK5/IJ ,-  M  '  KTOT  .  KMAX  r  KMAX2/BLK7/ ALF 
ll^u.O!-£CiA/,;I..K,'5/I.OTL(-50) 

Uj.;:c;,-4Sio  1  / '  3 f  c ( 3hi ) .  -jv 
1X(3)  .^(co)  f  r  (.■•:.K  -.AXi;)  »i.:. 

t-.vTLi.iiAL    '.jVAo._0 

KPior  -  (KT.n/2)+i 
Oj  o  I  rKi'.^x.  it!n;-iax 
A.:-t''ii:  -■  yo. 

;-.r3,.  iSi;    Ii.-ITIAL    COMOlTIOt, 
ii-     Cl.C'TtU.l)  )-  \,.illE    (6, 

I- jK>.Ar   (/'  t-:v(  •  f  ic^f ' )   - 


(  ITuvAX)  -P(3»3»3'51)  ?:<0KTli(35U  , 
L.PAK(7,7U1)  (YINr(&)  »0Yiirr(6) 


AR  I  A 
C  (  1  'I 


•  Z) 

LFAI' 
»70i.  )  . 


10 


(2i»::r->23,2'i  t2S,2 
-  AiiT'\G(r.v'(  i  I'cFn  ) 


i'^f.ALtuV  (ITLP)  ) 
Aj.,-;Ae(LV(  ITEK)  ) 


0,1 

Ai.FAR 

0  -l.&A 

G,j    To 

Ai>A!< 

AuFAi 

Go    To    1 

AlFAK    -    R£A1.  tEV(ITLR)  ) 

Rl    =    AlNAu ClV(ITEF) ) 

GO    To    1 

0,,£;GA    =    RlAL(LV(ntF')  ) 

R,:  -  Ai;''AG(e.V(rfi-:f;.) ) 

Go    To    ^ 

o.,t:GA  ■-  R;:/,i.(LV(nfii<) ) 
AlFai  "-  aa;iag(ev(itlr)  ) 

Go    To    1 
26    Rr    ::    i:FAl..  (r.y(  ITER)  ) 

ALFAi    =    AI.iAG(F.V  (  IIKR)  ) 

i    AuKA    :-    C;.,HLX(ALFAf;.ALFA 
Ir     (ALf'AP,     .LT.     0,0     .OU. 

2  GA'V;-:'\  -  CbryRT(ALFA^i'f^'>( 
SiG'-'A  =  CSo;;T(ALFA'fT;.:+( 
If  (!:'jTl-:(14)  )  LKllL'.  ih, 
F.,K;-:,>T  (  •  O'-.MV.A  -•  »it£i 
IF  (r<tAL(bA  ■i:A)  .LT.  U . 
Ir  (RtAL(brc:;A)  .LT.  C. 
IF  (;;01l126)  ./vi;C',  L'OTr: 
1  U'^P A! ;  C  2  ,  ,<  TOT  )  -C  ;•.£■_;.' ../A  L  F 

IF   (.mOT.   ,jotf(,;.o)    .c/;. 


5  AT  BOUNDARY  LAYER  EDGE 

130)    ii"ER»!-:y(rri£R) 

•  »r2E  1.6.8) 
6)  mXCV 


^2 


23 


."•■'j- 


<cO 


lijH 


Z(    L  I  i  tf-Pl    .,i   ) 

Zi;.f  if^PT'.  T) 

/I  3.  I  ;l'.moi') 

2(i»2»iNpror) 

Z^2f2/|^;''l".;^  ) 
ZC3»?.  fPiyl  ) 
Hi  /Of  K;>f,;  r ) 

^(3«  j»KFT.^I  ) 


(  1  .  I.W  0  .  !■ 

( 1 , ;)  /  u ,  u 

(  i)  .  U  f  1  ,  u 
-•  (  f  r  .  0  » I  . 

-•;■■.  0  >  I , 
„Lr  A*  >:'.<■ 
i  fi  .  0  .  0  .  ij 
(O.i'fil.-; 
(1 ,  0  »  0  .  L' 


1) 

0,;EGA    ,LT.     0.0     .OH.    RE    .LT,     C 
.0/1.0)  *1',F  *■  (  ALFAi  8LPAR  (2  /  KTOT  ) 
,Cp1  .  )  iRr:H';U,;'i- MALFA-n)LPAR(2fK 
liFl)     GAF'-'-'wSIGi-'A 
o  .,=. »  lOX;  'Si  OKA    -'  t2Fl6.fJ/) 
0)     G.'FiiA    -    "GAM-IA 
n )    'jIC-'A"-    --SIG;-;!\ 

; 3.J);  00  -   r>: ;iO(K roD-c (7.KT01 

A)/C(l/KroT) 

.I'jwT.   tj0ii-.'(3t;i )   c:j  -   (Or 0/0,0) 
1 
) 

)  -• :,  I  ;  =  ■  :A  (-Ol/  {  S I  GF,A -i  *.'i" ALf  A*  1 2  ) 
OJl-AIJ-^ 

■J) +ga:.!;.;a 

jU/  (SIGF.A  *  V2-ALF  A**2  ) 

) 

) 

) 


0)  R 
OMFG 
TOT) 


FTURM    9 
A)) 

-o;;t:GA) ) 


))  M 
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b4  DO    2    o'-i.  ril 

bb  /.iJfrii'^PJji)     =     (GA.vMAkZC  Jj  JL .  KPTOT  )  ~  (  0  .  .  I  .  )  *  7.  (  J  '  2  »  KPTOT  )  )  >  (ALFA-<- 

b6  IGai'-'  A)/Kr.~(0.  (1.  )  +;j|/<- (  4  .  -^aLFA/o. -!  G,\'.'rA*  (  ALt^  A+bA'-'MA  )  /  (  nLI-rtSlG'-.A)  ) 

bV  LZ;J(J>;^PTO"r  )/!■!£ 

bo  Z{.Ji':,y~'.-\.j-()     -    -C  (0. '1.  )  ^aLFA'GA'  ■1'**Z(  Jr  ).»!<PTCT)  +  (  ALFA  VGA!  i;A)  + 

59  .U(o.^-.KPTOT)  )  ^AL^'.»  (ALFAi-oA'''i'A)  .-.OL. '?  (  J  >  o  m<P  TOT  )  /  (  AL.F  A  +  S  I'isA  ) 
bO  0    /.  (w»o.rM"'Tv,l  )     =    -SIGI''.A*/:(  J;  J.-K;  TOT) 

61  Ir     (i.OTKCi'.  )  )    1,(.;ITL:    (iflOS)     C  (>:(J- I.KFTOT  )  r  I-)..;.;)  ?  J-l,;j) 

6^  li:b    F.jf<;-:Ar     (•     ImUIAL    VAl-UuS    JtlFOKt:    ooL.     II;TF..^'/{l^t..lO.'t)  ) 

60  1 1    -    1 

64  C.    IhT;  G.aaT;-    i.'ISTiJKt?AriCF    EOliS    FOR    AOF-U'ltu    EIGEf.'VAI-UEr  0i;TH0l!0i-;v,AL12It)G 

65  C    i-.Mf.K    ,.!LCt;b-)AKY    OK    '.v^•;£(<£    SPc.CIFIED 
6o  Ijo    t    KPp=ifKPTOT 

67  Kl--    =    KPTar+l-i<PP 

66  IF  (KP  .tlO.  KHTOT)  GO  TO  15 
69  Do  10  Jr-l.ii 

7U  K    =    2+KP+l 

71  D.J    9    !"1>M 

7a  y  Y.:(,T(i)   -  ,i(j»ipKp+i) 

7i  CmLL    LiIi-FcO{STAIit:u>DYIKT»  YIMTfil(K) 

7t  Dj    1m    1  =  1  f,M 

75  10    Z(J.I>*^P)    -    YINT(I) 

76  I'j    CaLL    OiUilCKCZflTLK/KPrArJGLEdlfPiKORTI-UAGLMIN) 

77  4    COuTli.UE 

76  IK     (t;0TE(17))    VvP.ITE    (6(117)     AGLMIIi 

79  H7    FOr^fsAT     {'    iMlNIf-^UM    ANGLE    BETV.'Etri    SOU^J    VECTORS    ENCOUNTEf^ED    III    ORTH, 

60  lP,iutE'->J    =•  (F10.6) 

61  IF  (AGLMI-i  .LT.  1.0)  RETU!?N  9 
02  5  IIEKu  -  U  -  1 

0^  IF  (r.0TE(17))  WRITE  (6.107)  IIENP 

8'(  iu7  Format  (•  total  ho   of  cian  reod  -'»i't) 

85  (,    EVALU,uE    bATt:-,FACTlOil    OF    B.C.    AT    Tl  iE    ivaLL 

66  Bi    =    i.LPAaCMl) 

87  02    =    (O.Gf(J.O) 

68  B3    =    (l.Of 0.0) 

69  X ( 1 )     =     ( 1  =  0 . 0 . 0 ) 

9J  DrTiO:    "     (  .,.' t  Z  (  J  r  &  ,  1  ) +33  *2  (  3  r  3  i  1  )  )  I  (  2  (  2  »  2  »  1 ) +B1  i^Z  (  2  ,  i  ,  1  )  ) - 

91  1(  jf.*Z(-'o<  1) +03  +  2(2.3/1)  )-f^(2(3.2f  J  )+U  1*2(3.1. 1)  ) 

92  Aid)    =■    X(^)  f  (  (b2'^2(1.6.1)+u3  +  2(1..3.  1)  )-5  (2(3.2.  l)+lUtZ(3r  1.  1)  )~ 
9  3  1  {^.':.  +  ..:(3,67  1  )+b3*Z(.3.3(  1)  )  ^  (2(1.?. ,  .1  )-(Bl*2(  1'  1  «  1)  )  )/0^:iJOM 

9'+  XC;)    -    .X  (  I  )  V  (  (1:2  +  2(2.6,  ].)+,i:H2(;:,;j,  i)  ) +  (2(1,1,  1  )+i!lt2(l.  1.  1)  )- 

9b  i(  ji*2(  J..'-.-  '.  )+rO'*''-  (  1  .o.l)  )♦  (2(2*2;  1  )  +E1  f-2  (2.  1.  1)  )  ) /UthlO.A 

9t,  D^j    In     X--:J  .    I 

97  lo    YU'l)     ■■     .(  1  )i<2(  1  .  .C  r  1  )+v  (2)  ■k2(2,  I  r  1  ) -i-X(,0  i2(3f  I  ,  1) 

V3  Ir     ('/./li:  (  i  v)  )     ..,..tlr.     (r..il2)      (  Y  (  I  .  .'. )  .  I- L  .  .1 ) 

99  li2    f  ...■.'■;. .\     (t     VC '1  );^' /12-:i0,") 

100  l,      (O  .':'.(  r('- .■  Ij  )     .lT.     l.Oc-2'j)    .^f-.TURN  ■ 

101  F(  ITL-J     -     I  (  1/  I  ) /Y (  i '  I  ) 

102  Ir     (.■:Al:i(;- (i  TEn)  !     -i-E.     .'"-.U-M    -A'-P.,     ITEP>     ,E(l.     1)     GO    Tq    7 

103  IF     (iA'f.:.  (.-:^i)  )  .,Pr:.-.  (6.1..■,^)  i  ■^■..,^  ,  F  (  ;  TT  ■( )  •  I  It; '1 ,  C  A«S  (  F  (  X  IT.  iO  ) 

lu'>  iCo  F-,-;i;;iaT  ('  r('»I--;.')  -  '  .  2i:io.'' . :  '-(r  '  Or.  CaUjC- (  '  '  i2)  '  )  ;  ;;'.L16.f,) 
1U5     C  !.;)Ti:-;,:':L  I;:-./!  .AL''!-.  O:'  AS.-.'JVc.j  ET.OtMVAl.iJE  TO  USE 

lot  I,"  '.  IT.-..:»-<)  .i'o,.} 

ij7  'i   caul  Ei.r.si(r  ,•  rn. ■'■/..■-.. I  V.  j'TE'oii.' ■.■p,.''.'..i,.Li;n) 
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lUfi  If    (r,uS(i,~!:EAL(Ev(n.;K)  )/.^!:alcev(  tTffc-n.) ) )    .GT.   t:PSL;j   .ci<. 

j.09  1A^S(  1  .~AI:  AGCiCv  (  I  ft!',)  )/Al.-lAG(£V  (  IT^Ril.)  )  )     .GT.    tPSLIs')    GO    TO    6 

lio  ir    (c;.i:  ^.(MITlR)  )    .Of.   cP'-iLri)   go   To  5 

lii  ir   ('.'Oir.(.:))   KdTurn 

1.U  U  J    .17    1-1  .3 

u^  17  xii)   =   ■;(  I  )/c,H<oS(>  (4..1) ) 

im  C-j    To    7 

Ub  6    CO,.TllJbE 

lltJ  IF     (liOT:.(2i))     i'.RITL     (f-fiCl)     EVCi;],Tfv)  .F('.:ITK) 

U7  r;^rc(£:X)    =    .FaLS£. 

ilii  ICi    F.jiv'-.A'r     (/'     ••AXIi-U;.     ALL.Gi/AulLF    riO    fTia;  AT  J  ONS    FOR    E  TGENVALUr;    EXCF  EPiflj 

ii.'-j  1'/'  iLSr  i::srif-:AT[:  of  e.v.  ;:•  »2ejh.3»  ' -attekpting  ro  satisfy  :hg'-:o  a 

iS^Cl  k.C.     v.'iT.i    /.u^-ITTaNCE    CF'»2E.l'Ub) 

-^21  IF   (iior::(.?2))   i.KnE   (otiii)    (ev(iter)  (F(iiEa)  ?  iter-i.  itr,^;ax) 

12^  lio    FjK.-iAf     (iJ.X'  'dV(  HlU)  '  .;'.5a.  'F(  1  fL'i.)  VClElG.S)  ) 

.'.i?3  Gj     fu    l.'>i'ii^»:^3f34».;/j»?,0)  »kEV 

lE'l-  ^1    AlFAR    -    i-'EALC£V(i;i  iR)  ) 

lc;b  0/:tGA    =    AiisAG(tV{f;iTI()  ) 

136  KcTURU 

12V  ;',2    AlFAR    =    »EAL(EVCrjI  it;)  ) 

i20  ALFA  I    =    Ai;,AG(cV(r;ITK)  ) 

i29  RETURN 

1^0  a    AlFAR    =    KEmL(EV((JITR)  ) 

131  Ht.    =    Al;iAG(EV(rJITR)  ) 

lo2  RETURN 

lii  :;H    0..,EGA    =    HEAL(EV(MITft)  ) 

li'+  RE    =    AlMAC,  (EVUaTK)  ) 

135  RETURN 

136  35    OiLGA    =    f;EAL(EV(UXTK)  ) 

137  A1.FAI    =    AI;-iA&(EV(filTR)  ) 
136  RETURN 

139  36    RE    =    HEAL((;V(I1ITR)  ) 

:!'+0  A1.FAI    =    A.t;;AG(EV(,,nU)  ) 

i'H  RETURN 

14.2  C    Ev/.LUhTE    £lGEiJFUflCTIO,-:S 

14.3  7    If     (NOTEdt))     V.uITE    (odOj)     ( Jt  X  ( J)  ,  J-l .  N) 

144  lu9    FuF-IAT     ('     VALUES    OF    COfijMl'iG    COEFF    AT    '.vALU:          X(';Il/')     ='»2tl5.3 

i't5  3/(3c.Xf 'xt  • »  Uf ' )   =',,:Ei').an 

I'+G  Call  eo:,fc;j{x. YtZi  i je:;.j-p»kortii) 

1^7  if     (,/jlE(,iJ))     ■.■.Mr-.     (o>U-''>)     (n.iU/iMTU  n)  '  n  =  l' I  J-LlfD) 

I'tB  114    r-Jn.-./T     ((I     :vC;:TH(  '  ,  :3t  '  )     -'fi'i)) 

l'*9  >'.,■(  HE     I  o.  10.2)     (isF»;.Li-A;^(  '. .  ,i  .Ki'-l )  .■  {  Y  (  I  (!',P)  »  [  =  1  '  ;■  )  ^K,P  =  1  ,,cr  TOT  ) 

150  1.j2    F  j,-;iAT     (/IjCX.  '  f  "  t,'j!)X*  •  »  = /c' <  »  '  KP  '  i  iVX  f  •  ET  A  '  f  6X  »  '  (R[iO  +  l 't  I  !)  '  »  .i'i  X  f  '  F  •  . 

ISL  lloX«  'Tl,  PuriATUitE  '  »  11  X;  '.--R^SSUr'E'  r  i;X  r  '  F  '  H  5X  p  '  TE'-T'ER  A  (  U'JE  '  /  (  I :»  i 

lb2  2F('.'i,12Eli),b)  ) 

lb3  iVi<lTEt'-5'  11))      (I'.p.rEi'A'.'C  r  L '"P  "1  '  .  t  <•.  ';;■'-.(  Y(  I  f  KP  )  )  .  I  -  1  >  "M  ,  ;<P=  I  »  KpTOT  ) 

IS'*  liC>    Fj!'.M,iT     C'     ,../.I..G     ihEiR    Cv:'.<!'ESFO.-;;:  j.;,G    A.';r'LITUJLS'/(  I'lM-y.'i  >6E.?0.e)  ) 

ib'J  RlIURN 

I'jG  e,,o 
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STAbEU 

1  SudKOuTln:    STAL-LQfY'OY^K! 

H  Coi-.rLtA     r?  jY»/.LFA>  .■.UciSSC 

i  L.jGICAL    i^dTE 

^  COMkOII    /tJLK4/'L;LHA>;iC/i:-LKb/rvrV.»KTor;Kf1AXrl<fIAX2/!'L.K7/ALFAK>  ALFaI  » 

!J  lf<c£»OvLv;A/,^LK6/r,OTi.;(30) 

6  DI^,E^'SIO,;    uY(M)  .YC:)  #^Li'Ai<(7r701 )  .-0(14.701) 

■7  AuFA    =    C  ii'LX(ALFA;;.ALFAn 
a  IF     (,.OT;„(co;  )     Cci'l    :=    ,1,0 

9  JK     (.i.'OT.    ,,OTL(<:^a))    COi.    =    0.0 

10  ir    (  .oif:(.::6) )   OEn  =  i.o 

11  If     (.;.C/T.    irOTE(,;6))    DEij    -    0.0 
l'^  .vuETjSC    -    iiLFAK(i.K)-OMiiGA/ALFA 

1-5  IF    (f-.OlE(-jO)  )    CY(i)=--{0,  f  X.)  MY{.-')KC(10m;)-C{7..  K)  )  =i  DEi-l «  ..LESSC*Y  (  J, 

14  1)  ) 

15  IF     (.ijOr,    ,jOTEC:50))     DY(1)     -    ~  (  O  ,  r  1  .  ) +r  (2) 

16  0H2.)    ~    Y(5) 

17  IF    (KOTE(iO))    DY(o)     ~    Y(6)*C0tl 

10  IF     (.l;UT.    i4OTE{30))    UY(3)    =    (  0  o  0  .•  J  .  0  ) 

19  IK     t(jOTr.(.iO)  )     DY(;)     =     {  0  .  r  1 .  )  t- (-ALFA**?.*V,'LF3riC*C  ( 1  »  K  )  *Y  (  I  )  +  (  ALFA/P. 

20  1E)*C(2'K)<  (--(0.  »1,  )>;(C  (o.K)-(ALFAiC(lri<)  )  *i=2  )  *  Y  (1 )  ■)  (  C  ( 'I » lO +C  (  ^  '  K  ) 

21  2+C(l'i<>^)  )  ^Y(2)-Y(L.)  I  (..LCSSC4DEr,MC(G»K)*Y(  j)+C(7»i-;)  ('Y(6)  ):CC(7'K) 

22  3*iJL(J-Cl;i',-s)  )  *i-,LPA,M  3i  K)  t.Y(3)  )  vCON)  ) 

23  IF     (.(lOr.    t.OTE(iO))    UYC^)    ::    (  0  .  p  1 .  )  i<  ( '-AlFA*  +  2  +  '.vLES5C*C  (1 » ;0  f  Y  ( 1 )  + 
^'^  lALFA/,-.E:Cl.i>K)  -f  (  (--Y(5)+Y{2)*C2.  ■'C(.5?K)+C(14fK)/2,  )  )/C(l.;^)-^- 
25  Z(u. »1. J+M.FA*+2*C(1»K)*Y(1) ) ) 

26  IF    (iloTrKoO))    DY(li)    ~    ALFa*RE  *■  (tjL.rAR  (3  >  K  ) +Y  {  1 )  i- (  0  ,  »  1 .  )  ^<  (  Y  (  ?.)  * 

27  lWi,ESSC  +  C(l.K)*Y(4)  )  )  *C  (  I .  K  ) /C  (2  >K)  ■!■(  ALr-A*C  (1  r  rJ  )*-+2*(  (  0  .  r  K  )  »C  ( 't' K 
2a  2  )  A Y  ( .1.  )  1  V  (  a )  ^  C  ( 1 U  .  K  )  -vv;i.EG.£C  + Y  (  3 )  iDEi4*C0N )  +C  ( 6 »  K  )  *Y  ( 5 )  +  ( C  (  9  *  K) 

29  Z>*)  {3)+C(L..i\)*Yl6)4;}LPAr'(3>K)  )  ^CON 

30  IF    (.iiOT.    i;01E(30))    r'Y(;>)    =    ALF-A^«r:*(Y(l)*ELPAR(3,K)  +  (0.  /  1  .  )  +  ( 

31  lY(2)+/.l-■;-:ji>C^C(l.K)-i■■YC'^)  )  )  i<C(l»K)  i  i  ^/C  C  2  r  K  ) +ALFA**2+ (  C  (1 »  K  )  **2*Y  (  2) 

32  2+  (  0  .  f  1  •  )  ^-  r  ( I )  )=  (  C  ( '•  '  ;<)  -C  (  !■+  r  K  )  /2  .  )  ;  !Y  (  S )  :^  ( C  ( C  »  K )  -C  ( 1'+  »  K  )  )  +  Y  ( 2  )  * 

33  3  ( C  ( 3  f  K  )  /  { ■:.  C  <+  »  K  )  -C  (  o  '  r. )  )  -C  (1 H  .  :•. )  )  .  C  (1 1!  ,  K  )  /2 . 

3'+  IF     (NOTEC^o)  )    (jY(<-,)     -     (ALFAt-KE:  CY  CD  "-hLPAKCbpK )  +  (().»  1.  ):-..LE3S.C*Yt  3 

35  1)  )*C(ilfK)+Y(3)*(  (hl.FA^CdfK)  )  *:*2  :-c  ( 12  »  \  )  ) +C  1 13  ,  l< )  :'Y  (  b)  )  iCOiJ 

36  IF    (.i.OT,    ,vOTE(30))    0Y(6)    =    (0.0.0.0) 

37  R£TUKri 
3a  E,,D 
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APPE>n)IX  E 

The  Coinpuirer  Prograiii  In  Operation  - 
Selection  of  Nuniarical  Parameters 

The  numerical  prop;rain  devalcped  to  solve  the  stability  pro- 
blem outlined  in  Chapters  II  and  III  is  written  in  FORTRAIJ  V  lan- 
guage and  evaluated  on  a  Univac  1103  digital  computer.   Just  as 
discussion  of  the  coraulete  problen  has  been  divided  into  fjo 
categories  (mean  and  disturbance) ,  so  too  is  the  computer  program: 

(1)  Solution  of  the  laean  flow  equations  in  real, 
double-precision  arithTuctic,  and  theii  cal- 
culation of  the  variable  coefficients  required 
for  solving  the  disturbance  equations; 

(2)  Solution  of  the  eigenvalue  problem  lis  complex, 
single-precision  arithmetic  and  subsequent 
manipulation  of  the  eigenf unctions  to  conipute 
Reynolds  stresses,  dissipation,  etc. 

Operation  of  all  or  part  of  these  two  eler.ients  is  regulated 
by  a  series  of  input  logical  variables. 

In  addition  to  the  parameters  (or  eigenvalues)  v/hich  expli- 
citly characterize  the  fonrulated  stability  problem  (a  ,  o.^  ^   w^, 
Re,  and  of  course  T  ,  T  ,  and  M) ,  there  are  also  adjustable  para- 
meters which  ijiplicitly  characterise  its  numerical  solution  (e.g., 
step  size,  orthonormalization  criterion,  specification  of  boundary 
layer  edge,  etc.).   Since  the  latter  govern  the   abiJity  to 
accurately  and  efficiently  determine  tlie  former,  discussDon  of  the 
program's  operation  will  be  developed  tiround  tlie  optim.um  selection 
of  these  numerical  paraiveters. 
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To  better  illustrate  the  difficulties  associated  with 
making  such  a  selection,  numerical  exanples  are  calculated  for  each 
parameter  for  the  same  typical  "eigenvalue"  x^fith  fixed: 

T  =  90°F,  T  =  60°F,  H  =  0,  a.,  =^  0,  and  Re  =  7000. 

\-7  CO  'I 

Table  E.l   demonstrates   the   sort  of   results   obtained   for   the   same 

fixed  parsr.eters   and  v/ith:      h,   =    .0175,    Q   --   88°,    rj  =   10,    and 

d        '  max 


n..  -  5.7925  (u.  =  .999020). 
&  o 


"r 

CO 

Coirments 

.123764 

3.72957 

including  fluctuating  temperature 
field  but  without  density 
fluctuations . 

.121905 

3.66240 

V'ithout  fluctuating  temperature 
field  but  with  all  r;ean  fluid 
properties  variable  in  dis- 
turbance equations. 

.122347 

3.68283 

Including  all  mean  and  fluctuating 
quantities,  equations  (2.37  -  2.40) 

.118809 

3.56218 

Without  fluctuating  temperature 
field  but  with  variable  mean 
viscosity  only  (Wazn;an.^  et.al 
equation  (1.2)  [3-6]). 

.119928 

3.60135 

Wazzan's  equation  using  property 
variation  described  by  Kaups  and 

Smith  [7]  (see  Figure  B.7). 

TABLE  E.l  Tlie  Influence  of  Included  Disturbance  Quantities  and 

Chosen  Property  Variation  on  an  l^ligenvalue  Calculation, 
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E.l  Boundary  Layer  Edpe,  ri    and  n^ 
1___J I;'ifJ5 !A 

Two  distinct  length  scales  must  be  specified  to  accurately 

solve  the  eigenvalue  problem:  n    ,  the  finite  value  of  the 

max 

independent  variable,  ri ,  at  which  the  mean  flow's  asyir:ptotic 
boundary  conditions  (eqns.  (2,21a))  are  assumed  to  be  approximately 
satisfied;  and  i-j.  (Oi-  <  ti    ),  arbitrarily  selected  boundary  layer 
edge  v/here  asymptotic  disturbance  flow  conditions  are  satisfied. 
For  computational  purposes,  there  is  no  reason  to  sup'^est  that  tlie 
two  need  be  the  same.   Riile  the  former  is  associated  with  obtain- 
ing accurate  mean  velocity  and  temperature  profiles,  the  latter 
determines  v;here  these  profiles  may  be  truncated  for  the  eigenvalue 
problem  without  appreciably  degrading  the  accuracy  of  its  solution. 

The  practical  capability'  for  determinine  the  prOj->er  n    , 
r  i  J  a      i   J     'max* 

however,  is  contingent  upon  utilization  of  an  efficient  means  for 

calculating  the  initially  unspecified,  raean-flov;,  wall  boundary 

conditions,  u'(0)  =  F"(0)  and  H.'(O).   The  effectiveness  of  the 

particular  method  used  in  this  investigation  [18]  for  finding  these 

conditions  is  demonstrated  by  the  example  in  Table  E.2  where 

T  =  90"?,  T  =  60"*?,  M  =  0,  h  =  .02,  p    =10. 
w       '   00  m        max 
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Iteration 

^'  STEP 

F"(0) 

H  '  (0) 

Guess 

1.0 

-1.0 

1 

2 

0.63890603 

-0.80589459 

2 

10 

0.45716795 

-0.71569664 

3 

10 

0.46637891 

-0.71328643 

A 

10 

0.46639972 

-0.71321361 

5 

10 

0.46639971 

-0.713213&1 

6 

10 

0.46639971 

-0.71321361 

TABLE  E.2   Converfrence  Hiscory  of  Unknown  Mean  Flow  Boundar>' 
Conditions,  r"(0)  and  ll'(O). 

With  0.55  seconds  required  for  the  first  iteration,  approximately 
2. 78  seconds/ iteration  thereafter,  and  six  iterations,  both  con- 
vergence to  values  of  F''(0)  and  H' (0)  with  better  than  eight-place 
accuracy  (for  the  specified  fluid  property  variation,  of  course),  and 
calculation  and  allocation  to  storage  of  all  tnean-flov;  coefficients 
needed  for  solution  of  the  disturbance  equation  can  be  accomplished 
in  less  than  twenty  seconds  (see  Table  E. 3). 

Final  selection  of  r\  is  determined  by  observing;  the  value 

ir.ax 

for  V7hich  an  increase  in  this  parameter  fails  to  n^ake  a  disccrnable 

chance  in  F"(0)  and  H'(0),   It  is  reasoned  that  values  of  n 

■"  ^  '  max 

greater  than  this  one  v.'ill  not  alter  the  mean  flow  profiles  (thus, 
the  variable  coefficients  for  the  disturbance  equations)  and  so 
must  represent  an  unnecessary  exnendjture  of  computational  effort 
and  tine. 
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The  results  cf  sucli  an  n    search  for  the  mean  vrall 

max 

temperature  and  velocity  gradients  considered  in  this  studj-  are 

presented  in  Table  E.3.   In  all  cases,  F"(0)  and  li'(O)  were  un~ 

chan,^ed  for  ri    >  10.   Thus,  notint?  that  the  time  required  to 
max  -  '       -' 

solve  the  mean  fJ.ow  equations  and  coitipute  disturbance  equation 

coefficients  varies  onlv  slightly  over  the  n    ranpe  considered, 

max     •  ' 

accuracy  rather  than  economy  governed  the  selection: 


n    H  10. 
max 

This  value  ensures  satisf actiofi  of  asymptotic  conditions  (eqn. 

(3.4))  to  order  (10   ).   Due  to  its  importance  in  calculating  the 

eigenvalue,  the  influence  of  p    on  p. ''^  is  also  displayed  in 
^  max     6  '  •' 

Table  E.3.   As  a  basis  of  comparison,  it  may  be  noted  that  the 

value  of  F"(0)  for  the  unheated,  flat  plate  case  agrees  well  with 

Howarth's  value  of  F"(0)  =  .33206  [9,  p.  129]. 

The  primary  considerations  for  selecting  p  V'ere  founded  in 

attainment  of  an  accurate  eip,enso].ution  and  reduction  of  runninj^ 

time  and  storage  requirements  (knov7n  a  priori  to  be  relatively 

large  when  eigenfunction  recovery  is  required  [20]).   Although  the 

time  restriction  could  be  relaxed  somewhat  if  the  information  were 

needed  badly  enough,  storage  limitations  (150,000  for  the  Univac 

8 

1108  computer)  could  not.   To  evaluate  the  effect  on  the  eigenvalues 

of  using  a  "truncated"  mean  flow  solution,  the  example  shown  in 

Table  E.4  vras  calculated  with:   T  =  90°F, 

w 
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n 

max 

F"(0)-u  (0) 

H'(0) 

n,=^   E 

(defined 
by  Eqn. 

Tim.e 
Required 
(sec) 

13.8464 

T  = 

V7 

60"? 

1.7152416 

.52435x10"^ 

6 

. 33255559 

7 

.33209439 

1.7203045 

.43557x10""^ 

15.3772 

8 

.33205907 

1.7207617 

.13211x10"^ 

16.7630 

9 

.33203738 

1.7207868 

.14727xl0"'^^ 

17.7298 

10 

.33205733 

1.7207876 

.6039  2xlO"-'-^ 

18.8970 

11 

.33205733 

1.7207876 

-20 
.91101x10 

20.9158 

T   = 

90°F 

6 

.^^6680267 

-.71341605 

1.5112973 

.19198x10"^ 

14.6256 

7 

.466<'!2702 

-.71322750 

1.5142917 

.12898x10"'' 

16.1574 

8 

.46640087 

-.71321419 

1.5145365 

.31727x10"-^° 

15.2932 

9 

.46639974 

-.71321362 

1.5145486 

.28695xlO~-'--^ 

17.0454 

10 

.46639971 

-.71321361 

1.5145491 

.95693xl0~-'-'' 

18.9308 

11 

.46639971 

-.71321360 

1.5145491 

.25249x10"^^ 

19.7002 

T 

=  150"^? 

-.78390563 

1.2098',  09 

.35870x10"^ 

14.5934 

6 

.74345286 

7 

.74320905 

-.78382177 

1.2111477 

.17598x10"^ 

14.5190 

8 

.74319434 

-.78381662 

1.2112361 

.31693xl0~^-'- 

15.8284 

9 

.74319379 

-.78381642 

1.2112399 

.21000x10"^'^ 

16.9238 

10 

.74319378 

-.78381642 

1.2112400 

.68194xl0""'-^ 

17.7392 

11 

.74319378 

-.78381642 

1.2112400 

.  17341x1 0"-'-''' 

19.3116 

TABLE  E.3  Influence  of  n  ,.  ,  on  computed  values  of  F"(0),  H'(0), 
and  Hr*. 

0 
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T   = 
w 

200°F 

6 

.96336141 

-.83923358 

1.0362211 

.11659x10"^ 

12 

3970 

7 

.96318594 

-.83918394 

1.0369269 

.47597x10"^ 

13 

5738 

8 

.96317728 

-.83918114 

1.0369723 

.713S8xlO~"'-^ 

14 

8100 

9 

.96317695 

-.83918104 

1.0369741 

.39403x10"-'^^ 

16. 

5640 

10 

.96317694 

-."83918104 

1.0369741 

.68671j:10~"'-^ 

18 

4704 

11 

.96317694 

-.83918104 

1.0369741 

.61905xlO~-'-^ 

18 

5032 

T  =  60"F,  M  =  0,  h 


m 


.02 


TABLE  E.3   Contini:ad 


T  =  60°F,  M  =  0,  1) 


.02,  f2  =  45",  Re.*  ==  7000,  a.,.  -  0  (lov/er 
0  i 


branch  without  density  fluctuations) 


^ 

"6    ' 

\   Steps 

#  Ortho- 
normali- 

Time/ 

«R 

T  rv6 

w  xlO 

zations 

Iteration 

7.00 

.999959 

351 

102 

11.1514 

.123762 

3.72956 

6.70 

.999903 

336 

102 

10.7356 

.123765 

3.72957 

5.80 

.999033 

291 

99 

9.5320 

.123817 

3.72989 

4.68 

.990055 

235 

96 

7.6570 

.124466 

3.73420 

3.14 

.900150 

158 

91 

5.1472 

.134641 

3.83808 

TABLE  K.4  Influence  of  the  Arbitrarily  Specified  Boundary  Layer 
Edge,  ru,  on  the  eip,envalue  computation. 

Note  that  consideration  of  additional  sejrnients  of  the 
boundary  layer  docs  not  appreciably  change  tlie  number  of  ortho- 
normalizations  required,  sup,Resting  that  most  of  thorn  are  required 
in  the  rc,':',ion  near  the  v^all.   It  does,  however,  require  more 
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Integration  steps  and  so  more  tir.-e  per  iteration.   Unfortunately, 
du&  to  storaf;e.  and  tin?e  liraitations,  it  is  more  difficult  to  as- 
cribe a  definite  value  to  rir  (also  the  car.c  for  the  disturbance 

0 

equation   step   size,   h  ,)    than   for  n  .    F.ven   usinp   the  ruaxiinuin 

■   d  TTsax 

allocated  storage  (h ,  =  .02  and  rir  =  7.0)}  variations  in  the 
eigenvalue  are  still  sli[^,htly  greater  than  their  indicated  sir. 
place  accuracy. 

Since  even  four  place  accuiracy  is  adequate  to  clearly 
delineate  between  calculations  with  and  without  temperature 
fluctuations  (see  example  in  section  3.4),  it  is  suffici.ent  to 
specify 

u^  =  .999 

0 

as  the  position  at  which  the  boundary  layer  may  be  truncated  for 

stability  calculations. 

Greater  accuracy  than  that  indicated  in  Table  E.4  is 

restricted  by  both  storage  and  the  interdependence  of  ri  and  h,  in 

the  program;  v;hen  h,  must  be  reduced,  so  too  must  n..  according  to 

the  relationyliip  h,  >  ni./350. 

d    0 
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E . 2  S tep  Size,  h 

Perhaps  no  other  parameter  more  critically  governs  the 
accuracy  and  efficiency  with  which  a  solution  is  nuniericall}'  cal- 
culated than  does  the  step  size.   If  it  were  chosen  too  large,  the 
fourth  order  Run^e-Kutta  and  Adaras-McuJ  ton  intep.ration  schenes 
(with  truncation  errors  proportional  to  h  )  would  he  basically 
inaccurate,  rapid  changes  of  the  dependent  variables  (such  as  in  the 
thermal  boundary  layer  for  large  Pr)  \/ould  be  obscured,  and  (most 
importantly  for  this  analysis)  between  consecutive  steps,  for  a 
given  aRe,  linear  independence  of  the  solution  vectors  may  be 
destroyed  beyond  accurate  restoration  by  the  orthonormalization 
procedure,  independent  of  the  fi  specified  (see  section  3.4.3).   If 
chosen  too  snail,  not  onjy  vjould  integration  time  and  storage  be 
prohibitively  large,  but  error  propagation,  enhanced  by  the  in- 
creased numbLT  of  mathematical  operations  required,  would  contami- 
nate the  solution  sooner  and  thus  require  more  frequent  ortho- 
normalization  within  a  given  interval  (see  section  3.3), 

Thus,  to  solve  the  disturbance  equations,  the  selected  h 
must  be  very  closely  related  to  the  specified  values  of  aRe  and  Q. 
Since  accurate  convergence  is  ensured  only  by  giving  the  ortho- 
normalization  procedure  a  chance  to  work  effectively,  as  aRe 
increases,  h  should  be  decreased  and  il   increased.   In  mapjiing  out 
a  stability  contour,  tiie  point  at  \.»lrich  convergence  difficxilty  can 
be  anticipated  is  preceded  by  a  limited  ranf';o  of  smaller  aRe  in 
which  the  number  of  item  Lions  reqijircd  for  convergence  is  greatly 
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Increased  and  more  importantly,  the  mininura  angle  encountered  be- 
tween solution  vectors  drastically  drops  to  less  than  1^ .   At  or 
beyond  this  point,  convergence  could  be  obtained  only  be  decreasing 
the  step  size.   It  is  interesting  to  note  that  similar  difficulties 
were  encountered  by  Hack  [24,  p.Z.33]  (in  a  schena  not  using 
orthoaonnallzation)  ,  but  at  nuch  smaller  values  of  a.Re. 

In  contrast,  for  the  numerical  values  calculated  using  the 
fourth  order  system,  of  Wazzan,  et.al  [3],  no  suclt  difficult}'  was 
found.   One  possible  explanation  can  be  advanced  by  considering 
the  examples  given  at  the  bef^inuing  of  section  3.4.   For  the  sixth 
order  system  (that  is,  the  problem  as  forraulated  in  this  study  but 
without  density  fluctuations) ,  the  mininium  angle  encountered  during 
the  boundary'  layer  integration  v7as  36.40'',  whereas  for  the  fourth 
order  system  (V^aszan's  equations)  it  was  78.21".   Since  the  region 
in  v.'hich  significant  thermal  fluctuations  exist  is  much  smaller 
than  that  of  the  velocity  fluctuations,  it  can  be  expected  that  the 
step  size  accordingly  must  be  smaller  for  the  sixth  order  system 
to  maintain  the  sar.e  degree  of  orthogonality  as  for  the  fourth. 

Since  the  Runge-Kutta  integration  scheme  requires  infor- 
mation at  lialf-step  intervals,  the  mean  eqxiations  are  integrated 
with  a  step  size,  h  ,  half  that  of  the  one  used  to  solve  the 
disturbance  equations,  h.,  (i.e.,  h   =  hj/2.).   This  permits  the 
program  to  calculate  the  required  variable,  mean-flow-dependent 
coefficients  for  the  latter  exactly  for  all  points. 
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An  example  illustrating  some  of   tlie  points   discitssed   above 
is   given   in   I'abJes   E.5   and  E.6   for: 


T   =  90°F,  T  =  60°F,  M  -  0, 

n    =  10,  "u,  =  .999,  f2  =  AS'",  Re,-  =  7000,  a^   ■■■^-   0  (  lower  branch 
max     '6  6  i 


without  density  fluctuations) 


h 

u'(0) 

u'(0) 

^6^^ 

.01 

.46639971 

-.71321361 

1.5145491 

.02 

.46639971 

-.71321361 

1.5145491 

.03 

.46639972 

-.71321358 

1.5145490 

.OA 

.46639973 

-.71321352 

1.5145489 

.05 

.46639977 

-.71321341 

1.5145-^189 

.06 

.46639982 

-.71321321 

1.5145487 

TABLE  E.5  Influence  of  Step  Size,  h  ,  on  the  Unspecified  Wall 


BT 


Temperature  and  Velocity  Gradients,  u'(0)  and  R  (0) 


\ 

//Steps 

Time/ 
Iteration(sec) 

j'/Orthonor- 
malizations 

«R 

to  xlO*^ 

CO 

.01675 

347 

11.1042 

103 

.123819 

3.72996 

.02 

291 

9.5320 

99 

.123817 

3.72989 

.03 

194 

6.4768 

94 

.123812 

3.72967 

.04 

146 

5.0016 

78 

,123797 

3.72917 

.05 

117 

4.0280 

67 

.123774 

3.72834 

.06 

98 

3.3430 

59 

.123741 

3.72718 

TABLE  E 

6  Influe 

nee 

of  the  Step 

Size 

,  h^j,  on 

the 

Coiiiputcd 

Eigen\ 

alue. 
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From  these  calculations  (aRe  =  866.6),  it  -.ppears  that  at  least 
four  place  eigenvalue  accuracy  is  obtainable  by  using  h  ,  <  .03. 
Note  that  this  choice  also  ensures  eight  place  accuracy  of  the 
mean  flow  solution  (for  the  specified  property  variation) .   For 
results  given  in  Chapter  IV>  a  step  size  of  h  =  .02  or  smaller  was 
used. 


E . 3  Orthouorroalization  Angle,  H. 

As  previously  indicated  (section  3.3),  the  orthonormali:?ation 
criterion  used  for  this  study  is  a  specified  mininurii  angle  per- 
mitted between  complex  solution  vectors.   (A  comparison  of  relative 
magnitudes  of  the  vectors  also  could  have  been  used  [39]).   If 
this  angle  v/ere  chosen  too  large  (e.g.,  90°  as  for  Wazzaii,  et.al. 
[3]),  aside  from  an  increased  integration  time,  the  extensive 
matrix  manipulation  required  for  the  orthonormalizatiou  procedure 
could  be  e:cpected  to  introduce  errors  in  the  solution  [39  ]. 
Alternatively,  if  it  were  chosen  too  small,  then  linear  independence 
of  solution  vectors  could  be  so  degraded  that  insufficient  resolu- 
tion would  remain  to  accurately  reestablish  an  orthonoriv.al  basis. 

Since  the  parai^ictric  group  aRe  appears  exponentially  in 
the  exact  solution  to  the  disturbance  equations  outside  the  boundary 
layer  (see  equations  (3.11  -  3.15))  or  inversely  as  a  coefficient 
mu].tiply in;z  the  hi^^est  order  derivative  ttn'm  (see  equations 
(2.2  -  1,A)),  it  can  be  anticipated  that  as  aRe  increases,  so  too 
should  the  rate  at  v.'hich  the  solution  vectors  linear  independence 
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degenerates.   It  follov;s,  then,  first  that  for  the  specified  n, 
the  number  of  orthoiiorrrializations  required  laust  also  increase. 
Secondly,  as  fi  is  increasing  to  accoi.imodate  the  larger  value  of 
aRe,  so  too  should  the  percentage  of  or thono realized  integration 
steps,  ranging;  from  zero  v;hen  Q  =  0°  to  one  hundred  when  fi  =  90°, 
These  observ-ations  are  substantiated  by  Figiire  E.l  and  Table  E.7. 

Care  must  be  exercised  to  ensure  t?!at  !^_  is  chosen  suffi- 
ciently large  that  the  resolution  of  solution  vectors  is  not  lost 
betv;een  successive  integration  steps  just  before  the  orthonor- 
malization  criterion  is  violated.   That  is,  if  the  smallest  angle 
between  solution  vectors  at  one  step  is  greater  than  .Q,  but  less 
than  n  at  the  next,  this  lunimura  angle  must  still  be  large  enough 
to  enable  reformation  of  the  orthonormal  basis.  \^en   degeneration 
of  orthogonality  is  suf ficientl}'  rapid  to  require  orthononnaliza- 
tion  at  consecutive  steps,  it  can  be  expected  that  the  minimum 
angle  between  solution  vectors  (and  thus  convergence  capability) 
v/ould  cease  to  be  a  function  of  Q.   and  become  one  of  the  step  size 
instead,  as  indicated  in  Figures  E.l  and  E.2. 
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Fif^ure  E.l   Influence  of  the  Orthonorma]  iz.-ition  Anr,ln,  J2,  on  the 
number  of  ortb.onornallzntions  required  per  iteration  and  on  the 
minimum  angle  encountered  in  the  larst  boundary  layer  inLep.ratioii. 
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9.0^ 


4J 


•a 
o 
u 

o 

o 
u 

w 


c 
< 

6 

E 


8.0   - 


w 
u 


5    7.0 


C3 

c 


6.0 


5.0  - 


A.O 


3.0   - 


2.0   - 


T     =    90°F 


GO°F 


u  -    .999 

0 

o 


Re 


a^=0        j 


lower  branch 


,02 


.04 


,06 


Fip.vire  E.2    Influence  of  the  step  sJze  on  the  ninimnm  angle 

encountered  betv;een  solution  vectors  on  the  last  eipenvalue 
Iteration  for  a  specified  il. 
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Nu 

nber  of 

n 

°R 

w  xlO^ 

CO 

Ortho 
tions 

loririalixa- 
required 

Minimum 
Angle 

5 

.179330 

3.70352 

4A 

1.3547 

10 

.179330 

3.70S52 

50 

1.5350 

20 

.179330 

3.70852 

74 

1.6S08 

30 

.179330 

3.70852 

79 

1.6808 

^0 

.179330 

3.70852 

83 

1.6808 

50 

.179330 

3.70852 

88 

1.6808 

60 

.179330 

3.70S52 

93 

1.6808 

70 

.179330 

3.70852 

100 

1.7733 

80 

.179330 

3.70852 

112 

11.3345 

88 

.179330 

3.70853 

140 

12.3368 

90 

.179330 

3.70853 

146 

not  computed 

where  Re  =  10  ,  a_  =  0,  b,  =  .04,  u  ■-=    .999  (upper  bra7ich,  r  t-   0) 

TABLE  E,7   Influence  of  Orthonormalizatiou  Angle,  Q,  on  Cov.iputed 
Eigenvalue. 

On  the  basis  of  these  calculations,  it  would  seem  that  al- 
though the  miniinuin  angle  between  solution  vectors  is  profoundly- 
affected  by  the  specified  orthonormalization  angle  in  some 
regions,  the  co'.7.puted  eigenvalue  is  not.   Further,  in  the 
intermediate  range  of  Q,    there  is  a  region  in  which  the  nur.lier 
of  required  orthonornalizations  is  significantly  less  than  for 
n  Oi  90°  and  v;here  this  nunber  varies  only  slightly  with  9.. 
Naturally,  such  a  reduction  translates  as  a  significant  savings 
in  required  progran  running  tire.   For  example,  for  the  Re  -■  7000 
(upper  branch)  calculation  sliown  in  Figure  E.l,  tlie  time  required 
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to  complete  eip,envalue  convergence  (in  seven  iterations)  is  re- 
duced fror!  36.9  Lo  34. A  scicor.ds  by  spacifying  Q   equal  50°  rather 
than  83°,  still  v;ith  alr.iost  six-place  accuracy.   Thus,  selection 
of  9.   apparently  can  be  based  primarily  on  computational  economy; 
that  is,  it  is  chosen  to  reduce  the  number  of  computer  manipulations 
and  program  storafre  requirements. 

Since  these  results  are  typical  of  others  found  during  the 
course  of  this  stiidy,  it  is  felt  that  the  selection 

n   =  45° 

represents  a  good  balance  betV7een  accuracy  and  efficiency.   To  note 
in  passing,  this  choice  is  also  compatible  with  the  results  of 
Gersting  and  Jankov/shi  [20]. 

E . 4   Ejpenvalue  Estimates,  a,  m,  and  Re 

As  indicated  in  section  3.3,  three  eigenvalue  estimates  are 
required  to  start  either  of  the  root-f indii'ig  schemes  used  in  this 
study  (i.e.,  that  of  Wazzan,  et.al.  or  Muller) .   To  trace  out  a 
contour  ir^  parameter  space,  tlie  first  eigenvalue  estimate,  ev  ,  is 
taken  from  the  lower  branch  of  the  stability  curves  of  Uazzan,  et.al, 
[4]  near  the  minimum  critical  Reynolds  number  (see  for  example 
Figure   4.4   ).   This  procedure  is  employed  for  two  reasons. 
First,  since  the  present  investigation  represents  a  departure  from 
the  Wazzan  results  due  to  the  inclusion  of  a  fluctuating  temperature 
field  and  consideration  of  all  fluid  property  variation,  it  wotild 
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seem  that  these  should  give  the  best  avnilable  initiitl  estimate. 

Second,  as  the  difficulty  of  solution  increases  v/lth  increasing 

aRe  (due  to  probleris  associated  with  degeneracy'  of  the  orthogonality 

of  solution  vectors,  etc.)?  ^-'^   eigenviilue  on  the  new  contour  is 

chosen  for  which  solution  should  be  easiest  (see  also  discussion  by 

Mack  [24,  p.  280]). 

The  second  and  third  eigenvalue  estimates  are  arbitrarily 

chosen  as  ev  -   1.01  ev  and  ev„  =  1.005  ev,„  +  1.02  ev   .   Uith 
Z  1        J  iK  11 

the  three  initial  eigenvalue  estimates  thus  specified,  the  three 
corresponding  boundary  layer  admittances,  Y  ,  can  be  calculated. 

The  plane  fitting  technique  of  VJas;.:an,  et.aJ.  is  used  first, 
conditionally  followed  by  that  of  Huller.   That  is  to  say,  only 
when  all  the  follovjing  criteria  are  violated  ^■/ill  the  latter  be 
utilized: 


ev . 

(1)     10"^  <  I— ^^-^^^ 1  <  10^ 

ev.^    ' 


(2)   from  the  denominator  of  a 

o 


0'  .Y  ,.  ,,  J^,  +  (Y  ,.  -,Y  .  y- 

;  oj  o(j-i)='^'  I    o(:i-2)  o]--'^'  I 


2 

■^  ^^o(i-l)\(j-2)^^ 


<  10-^-^ 
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(3)   iluller's  method  was  not  used  in  computlnp  the 

preceding  estimate.   (It  was  found  that  consecutive 
application  of  this  method  gave  poor  convergence.) 

In  general,  this  scheme  worked  v/ell,  ensuring  six  place 
accuracy  within  five  to  eight  iterations  (depending  on  the 
accuracy  of  the  initial  guess).  At  first,  to  ensure  uniform  con- 
vergence (after  the  first  three  arbitrary  estinates,  ev^  ->  ev, ) , 
an  estimated  "eigenvalue"  obtained  by  either  method  was  first 
averaged  witli  the  most  accurate  one  previously  calculated  (i.e., 
that  one  which  gave  the  smallest  value  of  |y  |)  before  reevaluating 
the  stability  equations.  However,  since  this  procedure  v;as  found 
to  increase  (by  as  much  as  50%)  the  number  of  iterations  required 
for  convergence,  it  was  modified  to  average  only  vjhen  the  condition 
of  uniform  convergence  was  violated  (i.e.,  when  the  condition 

Iy  ,.  ^,1  >  JY  .!  was  not  satisfied).   V.'hen  Muller's  method  v.'as  used 
'  o(i-l)  I    '  oj  ' 

alternately  as  discussed,  it  v;as  found  to  be  ver}'  effective,  usually 
refining  the  "eigenvalue"  estimate  so  as  to  reduce  |Y  |  by  several 
orders  of  magnitude.   Thus,  the  net  result  is  an  iteration  sclieme 
which  ultimately  but  sporntically  converges,  often  in  a  nonniontonic 
fashion. 

An  example  of  the  convergence  is  given  below  for  T^  =  90°!^, 
T  =  60°F,  M  =  0,  n  -  50",  h,  =  .OA,  u.  =  .999,  Re  =  7000,  a,  =  0, 
a„   =  .12233,  w  =  3.0821x10   (lov/er  branch  including  density 
fluctuations,  r  /  0) . 
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Iteration 

"r 

(J) 

Y 

0 

Guess  1 

.081212293 

.017239455 

.  432x10"" -'■ 

Guess  2 

.08202A413 

.017411849 

.342x10"^ 

Guess  3 

.08161835A 

.017584244 

.139 

1 

.081367993 

.017209441 

.176x10"-'- 

2 

.080728449 

.017008234 

.332x10"-^ 

3 

.080769324 

.017018352 

.323x10"^^ 

4 

.0S076913S 

.017018221 

.311x10"^ 

next  estimate 

.080769140 

.017013221 



TABLE  E.o   Conversence  History  of  an  2 


igenvalue  Search 


In  keeping  with  previous  discussion  of  obtainable  and  re- 
quired accuracy,  iteration  was  terininatcd  when  both 


Y 


\ 


and 


ev 


1  - 


JR 


ev 


(j  +  l)U 


and 


ev 


1 


iL 


ev 


(j+l)I 


<  10 


-5 


This  criterion  usually  ensured  "eigenvalue"  accuracy  to  six  places 
and  "zero"  wall  boundary  conditions  of  order  (10~^)  or  less. 
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